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ABSTRACT 

A  finite  element  formulation  was  used  to  solve  the  steady- 
state  two-dimensional  conduction  heat  transfer  equation  in 
the  condenser  wall  section  of  an  internally  finned  rotating 
heat  pipe.   A  FORTRAN  program  using  this  method  was  coupled 
with  the  COPES/CONMIN  program  for  sensitivity  analysis  of 
design  variables  and  for  automated  design  of  the  internal 
heat  pipe  geometry. 

With  water  as  the  working  fluid,  numerical  results  ob- 
tained for  copper  and  stainless  steel  heat  pipe  condenser 
sections  indicated  that  for  the  maximum  heat  transfer  rate, 
the  designer  should  machine  as  many  fins  as  the  condenser 
material  and  the  manufacturing  process  will  allow.   A  saw 
tooth  profile  is  preferable  to  spacing  between  fins. 
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I.   INTRODUCTION 

A.   THE  ROTATING  HEAT  PIPE 

The  rotating  heat  pipe  is  a  closed  container  designed 
to  transfer  a  large  amount  of  heat  in  rotating  machinery. 
Its  three  main  component  parts  are:   a  cylindrical  evapora- 
tor, a  truncated  cone  condenser,  and  a  working  fluid  as 
shown  in  Figure  1. 

At  rotation  above  the  critical  speed  of  a  rotating  heat 
pipe,  the  working  fluid  forms  an  annulus  in  the  evaporator, 
and  will  be  vaporized  by  heat  addition  to  it.   The  vapor 
flows  toward  the  condenser  as  a  result  of  a  pressure  differ- 
ence, transporting  the  latent  heat  of  vaporization  with  it. 
External  cooling  of  the  condenser  causes  the  vapor  to  con- 
dense on  the  inner  wall  and  release  its  latent  heat  of  eva- 
poration.  The  centrifugal  force  due  to  the  rotation  has  a 
component  acting  along  the  condenser  wall  that  will  act  to 
drive  the  condensate  back  to  the  evaporator  where  the  cycle 
is  repeated. 

In  a  conventional  heat  pipe,  the  force  driving  the  con- 
censate  back  to  the  evaporator  is  due  to  capillary  action, 
which  poses  a  limit  to  its  operation.   The  rotating  heat 
pipe  is  not  limited  by  capillary  action  and,  unlike  the 
thermosyphon  which  depends  on  gravity  to  cause  condensate 
return,  can  be  used  in  any  orientation  [1]. 


14 


B.   OPERATING  LIMITS  OF  A  ROTATING  HEAT  PIPE 

The  first  theoretical  investigation  of  the  rotating 
heat  pipe  conducted  at  the  Naval  Postgraduate  School  was  per- 
formed by  Ballback  [2]  in  1969.   He  studied  the  limitations 
of  performance  imposed  on  the  rotating  heat  pipe  due  to 
various  fluid  dynamic  mechanisms.   Using  existing  theory  and 
experimental  correlations,  he  was  able  to  estimate  the  sonic 
limit,  boiling  limit,  entrainment  limit,  and  the  condensing 
limit  of  performance. 

1 .   The  Sonic  Limit 

When  increasing  the  heat  flux  in  a  rotating  heat 
pipe,  it  is  possible  to  reach  a  limiting  flow  rate  of  the 
vapor  brought  on  by  a  choked  flow  condition  in  the  pipe. 
This  condition  imposes  a  limiting  value  on  the  amount  of 
energy  the  vapor  can  transport,  thus  reducing  the  effective- 
ness of  the  heat  pipe.   The  limiting  heat  transfer  rate 
becomes 


Q^  =  p   U  Ahf  (1) 

xt    v   v    fg 

and  the  vapor  velocity  is  considered  to  be  sonic, 


Uv  =  c  =  /  gQkRT  (2) 


where 


U   =  velocity  of  the  vapor  in  ft/sec,  and 

A   =  cross  sectional  area  for  the  vapor  flow  in  ft 


2 
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c    =  sonic  velocity  in  ft/sec 

2 

gn   =  gravitational  constant,  32.1739  f t-lbm/lbf-sec 

k    =  ratio  of  specific  heats 

R    =  gas  constant  in  ft-lbf/lbm  ?.R»  and 

T    =  absolute  temperature  in  ?.R. 

2 .   Boiling  Limit 

Kutateladze  [3]  postulated  that  the  transition  from 
nucleate  to  film  boiling  is  totally  a  hydrodynamic  process . 
He  determined  a  theoretical  formula  for  predicting  the  burn- 
out flux 

Qt  =  K  /p;  Ab  hfg{a  g(pf  -  PV))1A  (3) 

where 

K    =  constant  value 

3 
p    =  density  of  the  vapor  in  lbm/ft 

2 
A,   =  heat  transfer  area  m  the  boiler  in  ft 

h^   =  latent  heat  of  vaporization  in  Btu/lbm 

TO" 
■*■  o 

a    =  surface  tension  in  Ibf/ft 

2 
g    =  acceleration  of  gravity  m  ft/hr 

3 
p.p   =  density  of  fluid  in  lbm/ft 

3 

p    =  density  of  vapor  in  lbm/ft 

The  experimental  data  obtained  by  Kutateladze  suggested  a 
value  for  K  in  the  range  of  0% 13  to  0.19. 
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3 .   Entrainment  Limit 

The  flooding  constraint  in  a  wickless  heat  pipe  was 
examined  by  Sakhuja  [4]  who  developed  the  correlation 


A   C   h^  /  gDCpT-pTp 
=  _x fg_* f   v   v  (u) 

{1  +  Cpv/pf)    } 


where 


Q     =  heat  transfer  rate  in  Btu/hr 

2 
A     =  flow  area  m  ft 
x 

C     =  dimensionless  constant,  0.725  for  tube  with  sharp 
edged  flange 

h^    =  latent  heat  of  vaporization  in  Btu/lbm 

2 
g     =  acceleration  due  to  gravity  in  ft/hr" 

D     =  inside  diameter  of  heat  pipe  in  ft 

3 
pf        -   density  of  the  fluid  in  Ibm/ft 

3 
p     =  density  of  the  vapor  in  lbm/ft 


4  .   Condensing  Limit 

Ballback  [2]  determined  the  condensation  solution 
for  a  rotating  heat  pipe  by  modeling  the  condenser  section 
of  a  rotating  heat  pipe  as  a  rotating  truncated  cone.   He 
developed  the  following  expression  for  the  condensation 
limit : 

9  k.  p_  oj2  h.  {T  -T  }3  4  4    |  £ 

A     -     r2     ff  fgS    W,      ,ro        T    '     . T*  „   3,4 

Qt    =    tt{j  2^ }       {LR0    +    Lsin    <j>]    -R      } 


p-    sm   4> 


(5) 
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where 

Q     =  total  heat  transfer  rate  in  Btu/hr 

kf    =  thermal  conductivity  of  the  condensate  film 
1  in  Btu/hr-ft-°F 

3 
P-    =  density  of  fluid  in  lbm/ft 


=  angular  velocity  in  1/hr 


hr.  =  latent  heat  of  vaporization  in  Btu/lbm 

T  =  saturation  temperature  in  °F 

T  =  inside  wall  temperature  in  ?'F 
w  r 

up  =  viscosity  of  fluid  in  lbm/ft-hr 


<f> 


-   half  cone  angle  in  degrees 


R     =  minimum  wall  radius  in  ft 
o 

L     =  length  along  the  wall  of  the  condenser  in  ft 
u^    =  viscosity  of  the  fluid  in  lbm/ft-sec 

The  condensing  limit  equation  is  a  function  of  the  geometry 
and  speed  of  the  rotating  heat  pipe,  and  the  physical  pro- 
perties of  the  working  fluid. 

Tantrakul  [5]  calculated  these  limitations  for  a  heat 
pipe  with  specific  physical  characteristics  as  shown  in 
Table  I,  with  the  results  shown  in  Figure  2. 

TABLE  I 
Specification  of  a  Typical  Rotating  Heat  Pipe 


Length 

14.000 

inches 

Minimum  diameter 

2.000 

inches 

Wall  thickness 

0.125 

inches 

Internal  half  angle 

1.000 

degree 

Rotating  speed 

2700 

RPM 

Obviously  from  the  results  in  Figure  2,  the  condensing  limit 
is  the  predominant  limitation  for  the  amount  of  heat  that  can 
be  transferred  from  the  heat  pipe.   However,  the  other 
limitations  may  become  important  as  the  heat  pipe  geometry  and 
operating  conditions  are  varied. 

In  order  to  augment  the  heat  transfer  capacity  of  the 
heat  pipe,  recent  efforts  have  been  aimed  at  raising  the  con- 
densing limit  line  which  may  be  accomplished  by: 

a.  a  high  value  of  cone  angle,  to  increase  the  centrifugal 
driving  force, 

b.  some  type  of  promoter  of  dropwise  condensation  to  in- 
crease the  value  of  the  inside  heat  transfer  coeffi- 
cient, h,  or 

c.  use  of  an  internally  finned  condenser  to  increase 
the  inner  wall  surface  area  and  the  value  of  h,  since 
the  presence  of  a  fin  will  decrease  the  effective 
condensate  film  thickness. 

A  high  value  of  cone  angle  means  a  departure  from  the 
cylinder  or  shaft  shape.   Since  the  principal  known  appli- 
cation for  the  rotating  heat  pipe  is  in  the  cooling  of  rota- 
ting machinery  this  approach  was  not  pursued.   Although 
effective  promoters  of  dropwise  condensation  exist,  none  as 
yet  can  be  considered  to  be  permanent,  and  this  approach  was 
likewise  ruled  out.   The  remaining  alternative,  using  inter- 
nal fins  in  the  condenser  section  to  raise  the  condensing 
limit,  was  seen  as  the  best  choice. 
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C.   ANALYSIS  OF  THE  INTERNALLY  FINNED  ROTATING  HEAT  PIPE 

Pursuing  the  addition  of  internal  fins  as  a  way  to  raise 
the  condensing  limit,  Schafer  [6]  developed  an  analytical 
model  for  a  heat  pipe  with  a  triangular  fin  profile  as  shown 
in  Figure  3.   He  assumed  one-dimensional  heat  conduction 
through  the  wall  and  fin.   Corley  [7]  for  this  same  case 
developed  a  two-dimensional  heat  conduction  model  using  a 
Finite  Element  Method,  and  also  assumed  a  parabolic  tempera- 
ture distribution  along  the  fin  surface.   His  results  indi- 
cated a  significant  improvement  in  heat  transfer  performance 
of  about  75%  above  that  predicted  by  the  one-dimensional  model 
of  Schafer  [6].   However,  Corley  [7]  cautiously  noted  a 
probable  error  of  5  0%  existed  at  the  fin  apex,  and  conse- 
quently mentioned  that  there  may  be  a  total  heat  transfer 
error  of  as  high  as  15%.   Tantrakul  [5]  modified  Corley' s 
computer  program  by  increasing  the  number  of  finite  elements 
in  order  to  minimize  the  heat  transfer  error  at  the  apex  of 
the  fin.   His  results  with  this  modification  converged  with 
the  results  of  Corley.   Purnomo  [1]  developed  a  two-dimen- 
sional Finite  Element  Method  solution  using  a  linear  trian- 
gular finite  element  model  as  shown  in  Figure  4-.   Purnomo '  s 
[1]  Finite  Element  Method  program  also  worked  and  converged. 
Purnomo ' s  [1]  code,  when  made  to  approach  the  geometry  of  a 
smooth  tube,  did  not  agree  with  the  analytical  and  experi- 
mental data  obtained  by  Schafer  [6]  for  a  smooth  tube.   This 
cast  doubts  about  the  validity  of  Purnomo f s  code. 
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The  parametric  studies  conducted  using  Purnomo ' s  code 
gave  no  clear  indication  of  the  best  condenser  geometry  to 
maximize  heat  transfer.   Also,  his  code  was  tedious  to  use 
and  required  numerous  runs  to  obtain  data  since  it  was 
written  to  perform  only  one  analysis  at  a  time. 

Clearly  a  computer  program  that  could  make  numerous  runs 
with  minimal  data  input  and  could  also  automatically  find 
improved  designs  would  be  valuable. 

D.   THESIS  OBJECTIVES 

The  objectives  of  this  thesis  were  therefore: 

1.  To  modify  Purnomo ' s  [1]  computer  program  so  that  it 
is  compatible  with  the  COPES/COMNIN  program  [8]  and 
can  be  used  for  analysis  and  automated  design  of 
rotating  heat  pipes  (internally  finned  or  smooth). 

2.  To  compare  results  using  Purnomo ' s  code  with  analyti- 
cal results  for  a  smooth  tube  obtained  by  Schafer  [6] 
to  determine  if  and  where  an  error  exists. 

3.  To  use  the  sensitivity  analysis  capability  of  COPES/ 
CONMIN  in  conjunction  with  the  modified  program  to 
study  heat  transfer  in  an  internally  finned  rotating 
heat  pipe . 

4.  To  use  the  resulting  program  to  obtain  an  optimum 
design  for  an  internally  finned  rotating  heat  pipe 
to  obtain  experimental  data  to  compare  with  the 
analytical  results . 

5.  To  use  the  resulting  program  to  obtain  numerical 
results  in  place  of  data  obtained  from  expensive 
experimental  operations. 
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II.   NUMERICAL  OPTIMIZATION 

A.   BACKGROUND 

Most  design  processes  require  the  minimization  or  maxi- 
mization of  some  parameter  which  may  be  called  the  design 
objective.   For  the  design  to  be  acceptable,  it  must  satisfy 
a  variety  of  physical,  aesthetic,  economic  and,  on  occasion, 
political  limitations  which  are  referred  to  here  as  design 
constraints.   While  part  of  the  design  problem  may  not  be 
easily  quantified,  most  of  the  design  criteria  can  be 
described  in  numerical  terms. 

To  the  extent  that  the  problem  can  be  stated  in  numerical 
terms,  a  computer  program  can  be  written  to  perform  the 
necessary  calculations.   For  this  reason,  computer  analysis 
is  commonplace  in  most  engineering  organizations.   For 
example,  in  structural  design  the  configuration,  materials, 
and  loads  may  be  defined  and  a  finite  element  analysis  com- 
puter code  is  used  to  calculate  stresses,  deflections,  and 
other  response  quantities  of  interest.   If  any  of  these 
parameters  are  not  within  prescribed  bounds,  the  engineer 
may  change  the  structural  member  sizes  and  rerun  the  program. 
The  computer  code  therefore  provides  only  the  analysis  of  a 
proposed  design,  with  the  engineer  making  the  actual  design 
decisions.   This  approach  to  design,  which  may  be  called 
computer-aided  design,  is  commonly  used  today. 
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Another  common  use  of  analysis  codes  is  in  tradeoff 
studies.   For  example,  an  aircraft  trajectory  analysis  code 
may  be  run  repetitively  for  several  payloads ,  calculating 
the  aircraft  range,  to  determine  the  range-payload  sensitivity, 

A  logical  extension  to  computer-aided  design  is  fully 
automated  design,  where  the  computer  makes  the  actual  design 
decisions,  or  performs  trade-off  studies  with  a  minimum  of 
man-machine  interaction  [9]. 

B.   CONSTRAINED  FUNCTION  MINIMIZATION  (CONMIN) 

Vanderplaats  [10]  developed  an  optimization  program  CONMIN 
capable  of  optimizing  a  very  wide  class  of  engineering  pro- 
blems.  CONMIN  is  a  FORTRAN  program,  in  subroutine  form,  that 
optimizes  a  multi-variable  function  subject  to  a  set  of  in- 
equality constraints  based  on  Zoutendijk's  [11]  method  of 
feasible  directions  [12]. 

Three  basic  definitions  are  required  to  discuss  the 
use  of  CONMIN: 

Design  Variables  -  Those  parameters  which  the  optimization 
program  is  permitted  to  change  in  order  to  improve  the 
design.   Design  variables  appear  only  on  the  right  hand 
side  on  an  equation  and  are  continuous. 
Design  Constraints  -  Any  parameter  which  must  not  exceed 
specified  bounds  for  the  design  to  be  acceptable. 
Design  constraints  may  be  linear  or  nonlinear,  implicit 
or  explicit,  but  they  must  be  continuous  functions  of 
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the  design  variables.   Design  constraints  appear  only 
on  the  left  side  of  the  equations. 
Objective  Function  -  The  parameter  which  is  going  to  be 
minimized  or  maximized  during  the  optimization  process. 
The  objective  function  may  be  linear  or  nonlinear, 
implicit  or  explicit,  and  must  be  a  continuous  function 
of  the  design  variables.   The  objective  function  usually 
appears  on  the  left  side  of  an  equation,  but  it  may 
appear  on  the  right  side  if  it  is  also  a  design  variable. 
Design  constraints  and  objective  functions  are  usually  inter- 
changeable . 

C.  CONTROL  PROGRAM  FOR  ENGINEERING  SYNTHESIS  (COPES) 
Recall  that  the  optimization  program,  CONMIN,  was  written 

in  subroutine  form.   Vanderplaats  [8]  has  developed  a  main 
program  to  simplify  the  use  of  CONMIN  and  to  further  aid  in 
the  optimization  process.   The  user  must  supply  an  analysis 
subroutine  named  ANALIZ .   What  follows  are  programming  guide- 
lines to  ensure  compatability  with  COPES. 

D.  PROGRAMMING  GUIDELINES 

In  developing  any  computer  code  for  engineering  analysis, 
it  is  prudent  to  write  the  code  in  such  a  way  that  it  is 
easily  coupled  to  a  general  synthesis  program  such  as  COPES. 
Therefore,  a  general  programming  practice  is  outlined  here 
which  in  no  way  inhibits  the  use  of  the  computer  program  in 
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its  traditional  role  as  an  analytic  tool,  but  allows  for 
simple  adaption  to  COPES.   This  approach  is  considered  good 
programming  practice  and  provides  considerable  flexibility 
of  design  options.   Only  five  basic  rules  must  be  followed: 
I .   Write  the  code  in  subroutine  form  with  the  primary 
routine  called  as  SUBROUTINE  ANALIZ ( ICALC) .   The 
name  ANALIZ  is  compatable  with  the  COPES  program  and 
ICALC  is  a  calculation  control.   Note  that  subroutine 
ANALIZ  may  call  numerous  other  subroutines  as  required 
to  perform  the  necessary  calculations. 
II.   Segment  the  program  into  INPUT,  EXECUTION,  and  OUTPUT. 
The  calculation  control,  ICALC,  will  determine  the 
portion  of  the  analysis  code  to  be  executed.   ICALC=1; 
the  program  reads  all  data  required  to  perform  the 
analysis.   Also,  any  initialization  of  constants 
which  will  be  used  repetively  during  execution  is 
done  here.   This  initial  input  information  is  printed 
here  for  later  reference  and  for  program  debugging. 
ICALC=2;  the  program  performs  the  execution  phase  of 
the  analysis  task.   No  data  reading  or  printing  is 
done  here,  except  on  user-defined  scratch  disc.   Data 
may  be  printed  here  during  program  debugging,  in 
which  case  it  should  be  controlled  by  a  print  con- 
trol parameter  which  is  read  during  input.   In  this 
way,  this  print  may  be  turned  off  after  the  program 
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is  debugged,  but  may  be  used  again  during  future 
program  expansion  debugging.   The  reason  that  print- 
ing is  not  allowed  during  execution  is  that  when 
optimization  is  being  done,  the  code  will  be  called 
many  times  with  ICALC=2,  resulting  in  voluminous 
print.   ICALC=3;  the  results  of  the  anlaysis  are 
printed.   Also  the  essential  input  parameters  which 
may  have  been  changed  during  optimization  should  be 
printed  here  for  easy  reference.   In  summary,  when: 
ICALC  =  1  Read  input  data. 
ICALC  =  2  Execute  the  analysis. 
ICALC  =  3  Print  the  results. 
III.   Store  all  parameters  which  may  be  design  variables, 
objective  functions  or  constraints  in  a  single 
labeled  common  block  called  GLOBCM.   The  order  in 
which  they  are  stored  is  arbitrary.   A  listing  of  the 
COPES  program  should  be  checked  to  see  how  many  para- 
meters may  be  stored  in  GLOBCM  (the  dimension  of 
ARRAY) .   Initial  distribution  of  COPES  allows  for  1500 
parameters . 
IV.   During  execution  or  output,  no  parameters  which  are 
read  during  input  should  be  updated.   For  example, 
if  variable  X  is  initialized  during  input,  the 
execution  segment  must  not  update  X  such  as  X=X  + 
3.2.   Instead  a  new  variable,  Y=X  +  3.2  must  be 
defined . 
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V.   Write  all  programs  in  standard  language,  avoiding 
machine  dependent  capabilities  such  as  seven  letter 
FORTRAN  names  (CDC).   While  this  guideline  is  not 
essential  to  the  use  of  the  analysis  code  within  the 
COPES  program,  it  makes  the  analysis  code  much  more 
transportable  between  different  computer  systems ,  a 
capability  which  easily  justifies  a  slight  reduction 
in  efficiency  on  a  given  machine. 
Adherence  to  these  guidelines  not  only  leads  to  a  more 
readable  and  machine  independent  computer  code,  but  allows 
this  code  to  be  coupled  to  the  COPES  program  without  modifi- 
cation. 

Having  written  the  analysis  code,  it  may  be  executed 
either  with  a  simple  main  program  or  within  the  COPES  pro- 
gram to  perform  the  analysis.   To  insure  that  guideline  IV 
is  followed,  the  following  main  test  program  is  recommended. 
Note  that  this  program  calls  ANALIZ  twice  with  ICALC=2  and 
ICALC=3,  to  show  that  the  same  result  is  obtained  repeti- 
tively . 


C    MAIN  PROGRAM  TO  CHECK  SUBROUTINE  ANALIZ. 
C    READ,  EXECUTE ,  AND  PRINT 

DO  10  ICALC=1,3 
10   CALL  ANALIZ  ( ICALC ) 

C    EXECUTE  AND  PRINT  AGAIN  TO  BE  SURE  THE  RESULTS 
C    DO  NOT  CHANGE 

DO  20  ICALC=2,3 
2  0   CALL  ANALIZ  (ICALC) 

STOP 

END 
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III.   FINITE  ELEMENT  SOLUTION 

A.   REVIEW  OF  THE  PREVIOUS  ANALYSIS 

Schafer  [6]  studied  the  one-dimensional  model  heat  trans- 
fer solution  and  Corley  [7]  studied  the  two-dimensional 
model  for  an  internally  finned  rotating  heat  pipe.   Both 
used  the  same  assumptions  and  boundary  conditions  based  upon 
the  analysis  of  Ballback  [2],  which  are  similar  to  those  used 
in  the  Nusselt  analysis  of  film  condensation  on  a  flat  wall. 
The  more  important  of  those  assumptions  are: 

1.  steady  state  operation, 

2.  film  condensation,  as  opposed  to  dropwise  condensation, 

3 .  laminar  flow  of  the  condensate  film  along  both  the 
fin  and  the  trough, 

M- .   static  balance  of  forces  within  the  condensate, 

5.  one-dimensional  conduction  heat  transfer  through 
the  film  thickness  (no  convective  heat  transfer  in 
the  condensate  film) , 

6.  no  liquid  -  vapor  interfacial  shear  forces > 

7.  no  condensate  subcoolings 

8.  zero  heat  flux  boundary  conditions  on  both  sides  of 
the  wall  section  (symmetry  conditions),  as  shown 

in  Figure  5  > 

9.  saturation  temperature  at  the  fin  apex5 
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10.  zero  film  thickness  at  the  fin  apex,  and 

11.  negligible  curvature  of  the  condenser  wall. 
Purnomo  [1]  developed  a  two-dimensional  Finite  Element 

Method  solution  using  a  linear  triangular  finite  element 
model  as  shown  in  Figure  4.  Purnomo  modif  ied  Corley's  assump- 
tion that  the  fin  apex  was  at  the  saturation  temperature  and 
allowed  the  value  of  the  temperature  at  the  apex  to  float. 
He  assumed  a  parabolic  temperature  distribution  along  the 
fin  surface . 

Purnomo ' s  statement  of  the  problem  for  the  formulation 
of  the  Finite  Element  Method  as  shown  in  Figure  5  is 


2      2 

^-y  +  ^=0  (6) 

3x    3y 


with  the  boundary  conditions 


3T 
a)   along  boundary  1,     -k  ■*—  =  h,  (T  -  T   .) 
43        J  an    1       sat 


b)   along  boundary  3,     -k  ~  =  h„  (T  -  T  ) 

d  n     2       °° 


9  T 
c)   along  boundaries  2  and  4,    r —  =  0 

3  n 


A  detailed  description  of  the  numerical  formulation  is  pre- 
sented in  his  thesis . 
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Purnomo ' s  computer  program  consisted  of  a  main  program 
and  three  subroutines; 

a)  the  routine  "COORD"  used  to  define  positions  of 
system  coordinate  points, 

b)  the  routine  "FORMAF"  used  to  formulate  the 
Finite  Element  Method  equations,  and 

c)  the  routine  "BANDEC"  as  an  equation  solver  for  a 
symmetric  matrix  which  has  been  transformed  into 
banded  form. 

B.   THE  COMPUTER  PROGRAM  DEVELOPMENT 

Purnomo ' s  [1]  two-dimensional  finite  element  program 
is  the  basis  for  the  present  analysis  code.   The  first 
task  undertaken  in  the  development  of  this  thesis  was  to 
check  Purnomo ' s  code  for  validity.   Theoretical  results 
using  Purnomo ' s  code  for  the  case  of  zero  fins  were  compared 
to  theoretical  results  and  experimental  data  for  a  smooth 
tube  obtained  by  Schafer  [6].   Purnomo ' s  results  exceeded 
Schafer's  results  by  a  factor  of  approximately  two.   In 
studying  this  discrepancy,  an  error  was  discovered  in 
Purnomo f s  code.   In  encoding  formula  (11-11)  [1]  for  mass 
flow  rate  in  FORTRAN,  the'  sin  $  '  term  was  dropped.   The 
correct  form  of  the  equation  is  shown  below. 

2    2  2 

p   oj    (Rn    +   xsin   cj>)5  *(x)      sin  c|>  „ 

M        (x)    = [5  *(x)e+5  *(X)    tan   a] 

3    Uf 

(7) 
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The  effect  of  this  error  was  that  the  film  thickness  along 
the  length  of  the  condenser  remained  very  small.   Purnomo ' s 
original  code  consequently  indicated  a  much  higher  value  of 
heat  transfer  rate  than  the  correct  analysis  did.   In  a 
correct  analysis  code,  the  condensate  film  thickness  grows 
continuously  until  the  condensate  reaches  the  evaporator.   A 
plot  of  the  condensate  film  thickness  calculations  from 
Purnomo ' s  original  code  for  the  case  of  a  smooth  tube,  his 
corrected  code,  and  Schafer's  program  for  a  smooth  tube  are 
shown  in   Figure  6.   Purnomo ' s  corrected  code  for  the  case  of 
zero  fins  agreed  to  within  8  percent  of  the  results  for 
heat  transfer  rate  using  Schafer's  [6]  analysis  for  a  smooth 
tube . 

Since  no  experimental  data  existed  for  further  comparison 
of  the  finned  model,  the  next  task  undertaken  was  to  adapt 
the  analysis  code  to  permit  automated  design  and  sensitivity 
analysis  using  COPES/CONMIN .   Many  modifications  were  made, 
some  of  which  are  mentioned  here.   The  program  was  rewritten 
in  subroutine  form  and  segmented  into  input,  execution,  and 
output  sections  to  make  it  compatible  with  COPES.   Since 
COPES  was  written  to  use  single  precision  mathematics  and 
the  analysis  code  uses  double  precision  to  allow  for  possible 
ill  conditioning,  the  subroutine  ANALIZ  also  makes  the  trans- 
formation from  single  to  double  precision.   The  initial  value 
of  the  film  thickness  is  calculated  within  the  code  based  on 
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a  formulation  by  Sparrow  and  Gregg  [13].   Previously  a 
constant  value  based  on  the  same  analysis  was  used. 

The  modified  code  can  be  used  alone  for  analysis  of  a 
given  geometry  or  it  can  be  used  in  conjunction  with  COPES/ 
CONMIN  for  a  single  analysis,  sensitivity  analysis,  or 
automated  design  (optimization) . 

A  listing  of  the  revised  computer  program  is  included 
as  Appendix  B. 

C.   SENSITIVITY  ANALYSIS 

The  sensitivity  analysis  feature  of  COPES/CONMIN  was  used 
to  obtain  data  for  the  heat  transfer  rate  Q  as  a  function 
of  various  design  variables.   All  design  variables  were  held 
constant  except  the  one  of  interest  which  was  varied  as 
specified  in  block  Q  of  the  COPES  input  data. 

For  example,  to  obtain  the  heat  transfer  rate  Q  for  26 
different  values  of  rotational  speed,  the  following  data 
input  is  required: 


$  Block  P 
1 
15 

$  Block  Q 
16  ,  26 

3000.,  100.,  200.,  300.,  400.,  500.,  600.,  700., 
800.,  900.,  1000.,  2000.,  3000.,  4000.,  5000.,  6000., 
7000.,  8000.,  9000.,  10000.,  11000.,  12000.,  13000., 
14000.,  15000.,  3000. 
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For  each  plot  only  the  variable  on  the  ordinate  is 
changed.   All  other  design  variables  remain  constant.   The 
basic  design  used  for  sensitivity  plots  has  the  design 
variable  values  shown  below: 


TABLE  II 

BASIC  DESIGN  FOR  SENSITIVITY  ANALYSIS 

condenser  length 

cone  half  angle 

condenser  radius 

wall  thickness 

fin  height 

speed  of  rotation 

saturation  temperature 

ambient  temperature 

outside  heat  transfer 
coefficient 

fin  half  angle 

ratio  of  trough  width 
to  fin  base  width 

number  of  fins 


CLI 

r 

9.0  inches 

CANGL 

r 

1.0  DEGREES 

RBASEI 

= 

0.775  inches 

THICKI 

= 

0.03125  inches 

B 

= 

0.025  inches 

RPM 

= 

3000  RPM 

TSAT 

= 

100°  F 

TINF 

= 

60°  F 

HINF 

= 

5000  Btu/HR  FT 

FANGL 

= 

10 

ZOA 

= 

12.8 

ZFIN 

= 

40.0 

Note:   In  the  code,  ZFIN  appears  on  the  left  hand  side  of 
an  equation  and  is  therefore  by  strict  definition  not  a 
"design  variable" .   ZFIN  is  calculated  from  the  values  of 
CBASE,  EZERO,  and  EPSO. 

D.   DESIGN  OPTIMIZATION 

Counting  parameters  such  as  external  heat  transfer 
coefficient  HINF,  there  are  thirteen  possible  design 
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variables .   Nine  of  these  are  geometric  or  functional  para- 
meters such  as  wall  thickness,  fin  height,  and  speed  of 
rotation.   The  design  variables,  possible  constraint  func- 
tions, and  the  objective  function  appear  in  the  Global 
common  block  GL0B1  in  the  code  and  are  listed  below  by  fortran 
name  for  clarity. 

DESIGN  VARIABLES 
CLI 
CANGL 
RBASEI 
R2I 

THICKI 
BFIN 
TZ 
TSS 
TINF 
HINF 
FANGL 
ZOA 
RPM 

CONSTRAINT  FUNCTIONS 
ZFIN 
BOA 
DIFF 

OBJECTIVE  FUNCTION 
QTOT 

There  are  a  wide  variety  of  design  problems  that  can  be 
pursued  with  the  code.   For  example  one  might  wish  to 
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determine  the  smallest  length  condenser  section  and  the  best 
internal  geometry  for  a  specified  heat  transfer  rate- 
perhaps  to  cool  an  electric  motor. 

To  maximize  heat  transfer  rate  through  the  condenser 
wall  the  designer  adds  a  number  of  fins  to  increase  the 
inside  surface  area.   As  more  and  more  fins  are  added 
however,  the  cross-sectional  area  for  conduction  through 
each  fin  is  decreased.   Also,  the  film  thickness  of  the 
condensate  in  the  though  increases  and  in  fact  could  com- 
pletely cover  the  fins  and  substantially  reduce  heat  trans- 
fer through  the  fin.   So  there  should  exist  some  optimum 
combination  of  number  of  fins,  fin  height,  fin  half  angle, 
and  ratio  of  trough  width  to  fin  width  that  will  permit 
maximum  heat  transfer  rate . 

The  design  study  undertaken  was  to  determine  the  fin 
height,  fin  half  angle,  and  fin  spacing  which  would  yield 
the  maximum  heat  transfer  possible.   The  design  variables 
then  were  BFIN,  FANGL,  and  ZOA.   Other  potential  design 
variables  were  held  constant.   The  objective  function  to  be 
maximized  was  QTOT ,  heat  transfer  rate  out  of  the  condenser. 

For  comparison,  the  theoretical  upper  limit  on  heat 
transfer  was  calculated  based  on  an  external  surface  tempera- 
ture equal  to  the  working  fluid  saturation  temperature.   This 
assumes  that  there  is  no  thermal  resistance  across  the 
condensate  and  the  condenser  wall.   When  this  upper  limit 
was  used,  the  maximum  heat  transfer  rate  was  predicted  to  be 
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6  3,3  22  BTU/HR  using  the  following  formula: 


Q     =  h  2iTr  L  (Twall  -  T«)  (8) 

xmax 


where 

h    =  outside  convective  heat  transfer  coefficient 

(5000  BTU/HR-FT2-°F) 
r    =  average  outside  radius  of  condenser  wall  (0.0056  ft) 
L    =  condenser  length  (0.75  feet) 
Twall=temperature  of  the  outside  wail  (100°F) 
T°°   =  ambient  temperature  (60°F) 

Certain  constraints  were  placed  on  the  design  based  on 
engineering  judgement.   For  example,  the  number  of  fins  was 
not  allowed  to  exceed  40  0  and  the  minimum  fin  half  angle 
allowed  was  10  degrees .   These  values  were  based  on  structural 
and  manufacturing  considerations. 
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IV.       RESULTS 

A.  INTRODUCTION 

The  computer  code  was  used  in  conjunction  with  the  COPES/ 
CONMIN  program  for  sensitivity  analysis  and  design  optimiza- 
tion for  maximum  heat  transfer  rate.   Numerical  results  are 
discussed  below. 

B.  SENSITIVITY  ANALYSIS 

As  shown  in  Figures  7  and  8 ,  the  heat  transfer  rate  in- 
creases for  an  independent  increase  of  the  design  variable 
of  interest,  but  levels  off  at  the  theoretical  maximum  heat 
transfer  rate  for  the  given  overall  geometry.   The  heat  trans- 
fer rate  then  is  limited  by  the  external  resistance  which  has 
become  the  controlling  factor. 

In  Figure  9,  the  heat  transfer  rate  increases  in  a  simi- 
lar way  for  an  increase  in  external  heat  transfer  coefficient. 
Again,  the  rate  of  increase  appears  to  lessen  for  external 
heat  transfer  coefficients  above  10,000  BTU/HR-FT  -°F  due  to 
other  limiting  resistances.   In  Figure  10  the  heat  transfer 
rate  is  observed  to  increase  linearly  with  condenser  length. 
This  is  expected  since  Q  is  a  function  of  the  area  and  for  small 
values  of  $  the  area  varies  directly  with  length.   In  Figure  11, 
Q  is  seen  to  rise  in  a  non-linear  manner  for  an  increase  in 
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cone  half  angle,  <J> .   This  is  because  the  internal  heat  trans- 
fer coefficient,  h.  is  a  non-linear  function  of  cone  half 
angle. 

Purnomo  [1]  concluded  that  the  heat  transfer  rate  continu- 
ously increased  as  the  fin  half  angle  decreased  and  that  this 
was  largely  a  result  of  the  fact  that  the  number  of  fins  in- 
creased at  the  same  time.   He  also  stated  that  the  increase 
in  heat  transfer  was  only  slight  when  the  fin  half  angle  was 
less  than  11  degrees.   This  is  somewhat  misleading  in  that  in 
his  analysis  the  number  of  fins  was  being  changed  with  every 
change  in  fin  half  angle.   This  author  has  drawn  a  different 
conclusion.   Figure  12  shows  a  plot  of  heat  transfer  rate  vs. 
varying  fin  half  angle  all  for  a  condenser  with  UO  fins.   The 
heat  transfer  rate  as  a  function  of  fin  half  angle  rises 
sharply  from  1-11  degrees  and  continues  to  rise,  but  less 
steeply,  as  the  fin  half  angle  increases  beyond  11  degrees. 
The  maximum  fin  half  angle  possible  with  M-0  fins  is  67.5 
degrees . 

C.   CONSTRAINED  OPTIMIZATION 

In  the  design  problem  undertaken  to  determine  the  optimum 
internal  geometry  for  maximum  heat  transfer,  numerous  runs 
were  made  for  condensers  made  of  copper,  stainless  steel,  and 
a  ceramic  material.   These  materials  have  thermal  conduc- 
tivity values  of  231,9,  and  1.0  BTU/HR-FT-°F  respectively. 
It  was  expected  that  for  each  material  a  different  optimum 


38 


design  would  emerge.   The  results,  therefore,  were  unexpected 
From  different  starting  points  (original  designs) ,  some  out- 
side the  feasible  region,  and  for  external  heat  transfer 
coefficients  from  1000-50,000  BTU/HR-FT  .°F,  the  same  optimum 
design  for  maximum  heat  transfer  was  reached,  which  is  out- 
lined in  Table  III  below.   Each  material,  of  course,  had  a 
different  heat  transfer  rate  even  though  the  geometry  for 
maximum  heat  transfer  was  the  same. 

TABLE  III 

CONSTRAINED  OPTIMIZATION  RESULTS  FOR  COPPER 
STAINLESS  STEEL,  AND  CERAMIC  MATERIAL 

Optimum  Design  for  All  Materials 

fin  height  0.023  inches 

fin  half  angle  10.0  degrees 

number  of  fins  400.0 

ratio  of  trough  width  0  .  5 
to  fin  base  width 


MATERIAL 

Copper 

Stainless  Steel 
Ceramic  Material 


external 
(BTU/HR-FT2-0 

T) 

HEAT 

TRANSFER  RATE 
(BTU/HR) 

1000 
50,000 

13,822 
296,560 

1000 
50,000 

10,600 
38,650 

1000 
50,000 

3844 
5237 
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D.   UNCONSTRAINED  OPTIMIZATION 

The  results  of  the  constrained  optimization  runs  were  so 
unexpected  that  it  was  decided  to  remove  the  constraint  on 
the  number  of  fins  and  the  trough  width  to  fin  width  ratio 
and  repeat  the  optimization  runs.   If  another  identical  design 
was  again  reached  for  all  materials ,  the  code  would  have  to  be 
considered  in  error.   The  results  are  presented  in  Table  IV. 

The  optimum  design  for  maximum  heat  transfer  was  different 
for  each  material  and  each  heat  transfer  coefficient.   For 
each  material,  the  optimum  fin  height  b  increases  as  the  ex- 
ternal heat  transfer  coefficient  increases.   This  provides 
more  fin  surface  area  to  increase  heat  transfer  rate  as  the 
outside  heat  transfer  resistance  no  longer  controls.   Despite 
a  fifty-fold  change  in  heat  transfer  coefficient,  each 
material  maintained  essentially  a  constant  5  */b  ratio.   The 
higher  conductivity  materials  required  less  exposed  fin  sur- 
face than  the  lower  conductivity  materials.   Copper,  for 
instance,  could  be  57  percent  covered  by  the  trough  condensate 
while  the  ceramic  material  was  only  14-19  percent  covered  for 
maximum  heat  transfer. 

In  each  material  the  number  of  fins  decreased  for  an 
increase  in  external  heat  transfer  coefficient.   This  provides 
more  space  in  the  troughs  to  carry  the  increased  condensate 
which  results.   Otherwise,  fin  performance  would  be  degraded. 
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For  an  increase  in  external  heat  transfer  coefficient  the 
highest  percentage  change  in  heat  transfer  rate  occurs  in  the 
material  with  the  highest  thermal  conductivity.   In  the 
ceramic  material,  for  instance,  a  fifty-fold  increase  in  ex- 
ternal heat  transfer  coefficient  results  in  only  a  36  percent 
increase  in  heat  transfer  rate.   In  the  copper  condenser,  the 
heat  transfer  rate  is  increased  by  a  factor  of  21.   This  shows 
the  strong  dominance  of  the  wall  resistance  to  heat  transfer 
in  the  ceramic  material. 

For  all  three  materials,  the  trough  is  essentially  eli- 
minated for  the  optimum  design.   Why  the  trough  is  not 
entirely  eliminated  in  the  stainless  steel  condenser  section 
is  not  clear. 

Because  the  manufacture  of  a  very  large  number  of  fins  in 
a  small  diameter  condenser  is  not  practical,  it  is  desirable 
to  use  the  constrained  optimization  results  as  a  design  guide- 
line . 

A  comparison  of  the  unconstrained  and  constrained  design 
for  a  copper  condenser  section,  for  instance,  shows  that 
using  the  more  realistic  number  of  400  fins  and  its  proper 
fin  height  b  instead  of  204-0  fins  results  in  a  heat  transfer 
rate  difference  of  only  0.6  percent,   This  is  shown  clearly 
in  Table  V. 
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TABLE  V 


Comparison  of  Unconstrained  and  Constrained  Optimum  Design  for 
a  Copper  Condenser  Section  with  External9Heat  Transfer  Coeffi- 


cient  of  50000   BTU/HR-FT  -°F. 


Number 

b  (inches) 

of  Fins 

constrained 

0.023 

400 

unconstrained 

0.019 

2040 

%Dif f erence 
Q  (BTU/HR)      in  Q 

296,560        

29  8,33  0        0.6% 


The  author's  conclusion  is  that  within  the  realm  of  what 
can  now  be  manufactured,  the  same  basic  design  is  best  for 
all  materials  regardless  of  the  external  heat  transfer 
coefficient.   That  is,  the  designer  should  machine  as  many 
fins  as  the  condenser  material  and  the  manufacturing  process 
will  allow.   The  code  should  be  used  to  determine  the  fin 
height  to  avoid  degrading  the  fin  efficiency  by  too  high  a 
level  of  condensate  in  the  trough.   As  seen  in  Table  V,  it 
may  be  possible  to  lower  the  constraint  on  the  number  of  fins 
more  than  once  and  compare  the  resulting  heat  transfer  rate 
to  achieve  an  effective  but  less  expensive  design  to  manu- 
facture . 

E.   A  CAUTIONARY  NOTE 

The  reader  should  be  cautioned  that  when  an  analysis  code 
based  on  assumptions  made  for  certain  conditions  is  linked  to 
an  optimizer,  the  code  may  change  the  geometry  or  other 
conditions  such  that  the  original  assumptions  are  no  longer 
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valid.   Judicious  use  of  side  constraints  and  engineering 
common  sense  can  prevent  disaster.   Numerical  solutions  must 
not  be  blindly  accepted. 

It  should  be  noted  that  all  calculations  in  this  thesis 
were  made  with  water  as  the  working  fluid.   Sensitivity  data 
were  run  with  an  original  rotational  speed  of  3000  RPM  while 
design  optimization  was  done  for  a  rotational  speed  of  3600 
RPM.   Different 'combinations  of  operating  parameters,  working 
fluids,  and  condenser  section  materials  may  not  behave  in  a 
predictable  manner.   While  certain  trends  are  shown  for  the 
particular  condenser  section  studied  here,  it  would  be 
preferable  to  modify  the  code  to  analyze   the  particular 
problem  at  hand  than  to  extrapolate  these  findings  to  too 
broad  an  application. 

While  a  design  improvement  is  almost  certain  using  an 
optimization  routine  like  COPES/CONMIN ,  there  is  no  guarantee 
that  a  global  optimum  has  been  found.   Therefore,  for  a  given 
design  problem,  the  designer  should  start  the  optimization 
process  from  several  different  original  designs  and  compare 
the  results. 
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V.   CONCLUSIONS 

1.  The  computer  code  is  valid  and  can  be  used  for  single 
analysis,  sensitivity  analysis,  and  automated  design  of  an 
internally  finned  rotating  heat  pipe.   It  can  also  be  used  to 
generate  numerical  data  to  find  a  correlation  useful  for 
design. 

2.  For  zero  fins,  the  code  converges  to  the  results  for 
a  smooth  tube  obtained  by  Schafer  [6]. 

3 .  The  heat  transfer  rate  increases  for  an  independent 
increase  in  fin  half  angle,  rotational  speed,  and  number  of 
fins  but  levels  off  at  the  theoretical  maximum  heat  transfer 
rate  for  the  heat  pipe.   The  heat  transfer  rate  continues  to 
increase  for  an  increase  in  condenser  length  and  cone  half 
angle . 

4.  Within  the  realm  of  known  materials,  maximum  heat 
transfer  occurs  for  the  same  fin  geometry  regardless  of  the 
external  heat  transfer  coefficient.   That  is,  for  a  given 
condenser  radius,  one  should  machine  as  many  fins  as  possible 
with  the  condenser  material  used  to  maximize  heat  transfer. 
The  code  should  be  used  to  determine  the  correct  fin  height 
to  avoid  degradation  of  fin  performance  by  too  high  a  level 
of  condensate  in  the  trough. 
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5.   When  the  constraint  on  minimum  fin  half  angle  is  re- 
moved the  optimum  design  is  a  large  number  of  very  thin  fins. 
A  practical  design  then  might  be  a  series  of  long  rectangular 
fins  . 
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VI.   RECOMMENDATIONS 

Based  upon  the  calculated  results  of  this  thesis,  the 
following  recommendations  are  made. 

1.  Build  an  internally  finned  heat  pipe  to  obtain 
experimental  data  for  comparison  to  the  analytical 
predictions . 

2.  Generate  a  new  model  to  analyze  the  case  where  cone 
half  angle  is  zero.   This  would  greatly  decrease  the 
cost  of  manufacture,  and  permit  the  use  of  internally 
finned  heat  exchanger  tubing  which  is  commercially 
available . 

3.  Analyze  different  shaped  fins,  including  rectangular. 
M- .   Obtain  numerical  data  using  the  code  to  attempt  to 

find  a  correlation  for  the  heat  transfer  rate  as  a 
function  of  important  variables: 

Q  =  Q  (temperature,  condenser  geometry,  fluid 

properties,  condenser  material  properties, 
heat  transfer  coefficients,  etc.) 
5.   Modify  the  code  to  include  as  constraint  functions 
structural  failure  modes  such  as  rupture  due  to  the 
internal  spinning  mass,  buckling,  and  whirling. 
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Figure  1.   Schematic  Drawing  of  a  Rotating  Heat  Pipe 
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Figure    2.      Operating   Limits    of   a  Typical,    Water-Filled 
Rotating   Heat    Pipe 
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Figure  3 .   Internally  Finned  Condenser  Geometry 
Showing  Fins,  Troughs,  and  Lines  of 
Symmetry 
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Figure  5.   Differential  Equation  and  Boundary  Conditions 
Considered  in  the  Analysis  of  Purnomo  [1] 
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Figure  7.   Heat  Transfer  Rate  (Q)  vs.  Number  of  Fins 
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Figure  10.   Heat  Transfer  Rate  (Q)  vs.  Condenser 
Length 


57 


120- 


100- 


?. 


X 

on 
X 

D 

i— 

CO 


80- 


LU 
< 


60- 


Z 
< 

a: 


< 

LU 

X 


40- 


20- 


~T~ 

5 


— r- 
10 


15 


— r— 
20 


CONE  HALF  ANGLE  (DEGREES^ 


Figure  11.   Heat  Transfer  Rate  (Q)  vs.  Cone  Half  Angle 
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APPENDIX  A 
USER'S  MANUAL 

This  appendix  describes  the  data  cards  for  the  use  of 
the  computer  program. 

The  data  is  divided  into  the  COPES/CONMIN  section  and  the 
heat  pipe  analysis  program  section. 

The  COPES  data  is  segmented  into  "blocks"  for  convenience 
All  formats  are  alphanumeric  for  TITLE  and  END  cards,  F10  for 
real  data,  and  110  for  integer  data.  Comment  cards  may  be 
inserted  anywhere  in  the  data  deck  prior  to  the  END  card  and 
are  identified  by  a  dollar  sign  ($)  in  column  one.  The  COPES 
data  deck  must  terminate  with  an  end  card  containing  the  word 
'END'  in  columns  1-3. 

Information  is  included  in  this  appendix  pertaining  to 
data  needed  for  single  analysis,  sensitivity  analysis,  or 
optimization  using  COPES/CONMIN. 

UNFORMATTED  DATA  INPUT 

While  the  user's  sheet  defines  COPES  data  in  formatted 
fields  of  ten,  the  data  may  actually  be  read  in  a  simplified 
fashion  by  separating  data  by  commas  or  one  or  more  blanks. 
If  more  than  one  number  is  contained  on  an  unformatted  data 
card,  a  comma  must  appear  somewhere  on  the  card.   If  exponen- 
tial numbers  such  as  2.5+10  are  read  on  an  unformatted  card, 
there  must  be  no  embedded  blanks.   Unformatted  cards  may  be 
intermingled  with  formatted  cards .   Real  numbers  on  an 
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unformatted  card  must  have  a  decimal  point. 
EXAMPLES 

Unformatted  data-, 
5,7,1.3,1.0+20,-5.1 
5,7,1.3,1.0  +  20,  ,-5.1 
5  7  1.3  1.0  +  20,  ,  -5.1 
5       7   1.3,  1.0+20  0  -5.1 
Equivalent  formatted  data; 

col-*  10      20      30      40      50      60      70      80 
5       7      1.3   1.0+20     0     -5.1 
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DATA  BLOCK   A 


DESCRIPTION:   COPES  TITLE  CARD 


FORMAT:     20A4 


Title 


FIELD  CONTENTS 

1-8  Any  80  Character  title  may  be  given  on  this  card 
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DATA  BLOCK   B 

DESCRIPTION:   PROGRAM  CONTROL  PARAMETERS 
7110 


NCALC    NDV     NSV 


CONTENTS 


1  NCALC:   Calculation  control 

0-  Read  input  and  stop.   Data  of  blocks  A,B 
and  V  is  required.   Remaining  data  is 
optional . 

1-  One  cycle  through  program.   The  same  as 
executing  ANALIZ  stand-alone.   Data  of 
blocks  A,B  and  V  is  required,  remaining 
data  is  optional. 

2-  Optimization.   Data  of  Blocks  of  A-I  and  V 
is  required.   Remaining  data  is  optional. 

3-  Sensitivity  analysis.   Data  of  blocks, 
A,B,P,Q  and  V  is  required.   Remaining  data 
is  optional. 

2  NDV:   Number  of  independent  design  variables  in 

optimization. 

3  NSV:   Number  of  variables  on  which  sensitivity 

analysis  will  be  performed. 
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DATA  BLOCK   C 


Omit  if  NDV=0  in  block  B 


DESCRIPTION:   INTEGER  OPTIMIZATION  CONTROL  PARAMETERS 
FORMAT:   7110 


IPRINT   ITMAX  ICNDIR   NSCAL   ITRM   LINOBJ  NACMXI 


FIELD 


CONTENTS 


IPRINT-  Print  control  used  in  the  optimization 
program  CONMIN. 

0-  No  print  during  optimization. 

1-  Print  initial  and  final  optimization 
information . 

2-  Print  above  plus  objective  function  value 
and  design  variable  values  at  each  iteration. 

3-  Print  above  plus  constraint  values,  direction 
vector  and  move  parameter  at  each  iteration. 

M--  Print  above  plus  gradient  information. 

5-  Print  above  plus  each  proposed  design  vector, 
objective  function  and  constraint  values 
during  the  one  dimensional  search. 

ITMAX:   Maximum  number  of  optimization  inter- 
actions allowed.   Default=20. 

ICNDIR:   Conjugate  direction  restart  parameter. 
GT.  0-  Scale  design  variables  to  order  of 
magnitude  one  every  NSCAL  iterations.   LT .  0- 
Scale  design  variables  according  to  user  in- 
put scaling  values.   If  not  zero,  NDV  +  1 
is  recommended. 

ITRM:   Number  of  consecutive  iterations  which 
must  satisfy  relative  or  absolute  convergence 
criterion  before  optimization  process  is 
terminated.   Default=3. 

LINOBJ:   Linear  objective  function  identifier. 
If  the  optimization  objective  is  known  to  be 
a  linear  function  of  the  design  variables, 
set  LIN0BJ=1.   Default=Non-linear . 
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NACMX1 :   One  plus  the  maximum  number  of  active 
constraints  anticipated.   Def ault-NDV-2 . 
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DATA  BLOCK   D     Omit  if  NDV=0  in  block  B 

DESCRIPTION:   FLOATING  POINT  OPTIMIZATION  PROGRAM  PARAMETERS 

FORMAT:     7F10 


FDCH    FDCHM    CT     CTMIN    CTL   CTLMIN   THETA 


NOTE:   Two  cards  are  read  here. 

FIELD  CONTENTS 

1  FDCH:   Relative  change  in  design  variables  in 

calculating  finite  difference  gradients . 
Default=0 .01. 

2  FDCHM:   Minimum  absolute  step  in  finite  dif- 

ference gradient  calculations.   Default= 
0.001. 

3  CT :   Constraint  thickness  parameter.   Default1 

=-0.05. 

4  CTMIN:   Minimum  absolute  value  of  CT  consi- 

dered in  the  optimization  process. 
Default=0.004. 

5  CTL:   Constraint  thickness  parameter  for 

linear  constraints.   Default=-0 . 01 . 

6  CTLMIN:   Minimum  absolute  value  of  CTL  con- 

sidered in  the  optimization  process. 
Default=0 .001. 

7  THETA:   Mean  value  of  push-off  factor  in  the 

method  of  feasible  directions.   Default= 
1.0. 
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DATA  BLOCK   D     Omit  if  NDV=0  in  block  B 
FORMAT:     4F10 


DELFUN  DABFUN  ALPHAX  ABOBJ1 


FIELD 


CONTENTS 

DELFUN:   Minimum  relative  change  in  objective 
function  to  indicate  convergence  of  the 
optimization  process.   DEFAULT=0 . 001 . 

DABFUN:   Minimum  absolute  change  in  objective 
function  to  indicate  convergence  of  the 
optimization  process.   DEFAULT=0 . 001  times 
the  initial  objective  value. 

ALPHAX:   Maximum  fractional  change  in  any  design 
variable  for  first  estimate  of  the  step  in 
the  one-dimensional  search.   DEFAULT=0.1. 

AB0BJ1:   Expected  fractional  change  in  the 

objective  function  for  first  estimate  of  the 
step  in  the  one-dimensional  search. 
DEFAULT=0.1. 


REMARKS 

1)   The  DEFAULT  values  for  these  parameters  usually  work  well 
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DATA  BLOCK   E     Omit  if  NDV=0  in  block  B 

DESCRIPTION:   TOTAL  NUMBER  OF  DESIGN  VARIABLES,  DESIGN 
OBJECTIVE  IDENTIFICATION  AND  SIGN 

FORMAT:     2110, F10 


NDVTOT   IOBJ   SGNOPT 


FIELD  CONTENTS 

1  NDVTOT:   Total  number  of  variables  linked  to  the 

design  variables.   This  option  allows  two  or 
more  parameters  to  be  assigned  to  a  single 
design  variable.   The  value  of  each  para- 
meter is  the  value  of  the  design  variable 
times  a  multiplier,  which  may  be  different 
for  each  parameter.   DEFAULT=NDV. 

2  IOBJ:   Global  variable  location  associated  with 

the  objective  function  in  optimization 

3  SGNOPT:   Sign  used  to  identify  whether  function 

is  to  be  maximized  or  minimized.   +1.0  indi- 
cates maximization.   -1.0  indicates  minimi- 
zation.  If  SGNOPT  is  not  unity  in  magnitude, 
it  acts  as  a  multiplier  as  well,  to  scale  the 
magnitude  of  the  objective. 
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DATA  BLOCK   F    Omit  if  NDV=0  in  block  B 

DESCRIPTION:   DESIGN  VARIABLE  BOUNDS,  INITIAL  VALUES  AND 
SCALING  FACTORS 

FORMAT:     4F10 


VLB     VUB      X     SCAL 


NOTE:   READ  ONE  CARD  FOR  EACH  OF  THE  NDV  INDEPENDENT  DESIGN 
VARIABLES. 

FIELD  CONTENTS 

1  VLB:   Lower  bound  on  the  design  variable.   If 

VLB.LT.-1.0E+15,  no  lower  bound. 

2  VUB:   Upper  bound  on  the  design  variable.   If 

VUB.GT.10.E+15 ,  no  upper  bound. 

3  X:   Initial  value  of  the  design  variable.   If 

X  is  non-zero,  this  will  supercede  the  value 
initialized  by  the  user-supplied  subroutine 
ANALIZ. 
4-  SCAL:   Design  variable  scale  factor.   Not  used 

if  NSCAL.GE.O  in  BLOCK  C. 
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DATA  BLOCK   G     Omit  if  NDV=0  in  block  B 
DESCRIPTION:   DESIGN  VARIABLE  IDENTIFICATION 
FORMAT :   2110, F10 

1       2       3       4       5       6 


NDSGN   IDSGN   AMULT 


NOTE:   READ  ONE  CARD  FOR  EACH  OF  THE  NDVTOT  DESIGN  VARIABLES 
IN  THE  SAME  ORDER  AS  IN  BLOCK  F. 


FIELD  CONTENTS 

1  NDSGN:   Design  variable  number  associated  with 

this  variable. 

2  IDSGN:   Global  variable  number  associated  with 

this  variable. 

3  AMULT:   Constant  multiplier  on  this  variable. 

The  value  of  the  variable  will  be  the  value 
of  the  design  variable,  NDSGN,  times  AMULT. 
DEFAULT =1.0. 
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DATA  BLOCK   H     Omit  IF=NDV  Q  in  block  B 
DESCRIPTION:   NUMBER  OF  CONSTRAINED  PARAMETERS 


FORMAT 


110 


NCONS 


FIELD 


CONTENTS 

NCONS:   Number  of  constraint  sets  in  the 
optimization  problem. 


REMARKS : 

1)   If  two  or  more  adjacent  parameters  in  the  global  common 
block  have  the  same  limits  imposed,  these  are  part  of  the 
same  constraint  set. 
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DATA  BLOCK:   I     Omit  if  NDV=0  in  block  B,  or  MCONS=0  in 

block  H 

DESCRIPTION:   CONSTRAINT  IDENTIFICATION  AND  CONSTRAINT  BOUNDS 

FORMAT:   3110 


ICON    JCON    LCON 


NOTE:   READ  TWO  CARDS  FOR  EACH  OF  THE  NCONS  CONSTRAINT  SETS. 

FIELD  CONTENTS 

1  ICON:   First  global  number  corresponding  to  the 

constraint  set. 

2  JCON:   Last  global  number  corresponding  to  the 

constraint  set.   DEFAULT=ICON. 

3  LCON:   Linear  constraint  identifier  for  this 

constraint  set.   LC0N=1  indicates  linear 
constraints . 
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DATA  BLOCK:   I     (CONTINUED) 


FORMAT:   4F10 


BL     SCAL1    BU     SCAL2 


FIELD  CONTENTS 

1  BL:   Lower  bound  on  the  constrained  variables 

If  BL.LT.-1.0E+15 ,  no  lower  bound. 

2  SCAL1:   Normalization  factor  on  lower  bound. 

DEFAULT=MAX  of  ABS(BL),  0.1. 

3  BU:   Upper  bound  on  the  constrained  variables 

If  BU.GT.1.0E+15,  no  upper  bound. 
4-  SCAL2:   Normalization  factor  on  upper  bound. 

DEFAULT=MAX  of  ABS(BU),  0.1. 


REMARKS : 

1)  The  normalization  factor  should  usually  be  defaulted. 

2)  The  constraint  functions  sent  to  CONMIN  are  of  the  form: 

(BL- VALUE )/SCALl  . LE .  0.0  and  (VALUE  -  BU)/SCAL2  . LE . 
0.0. 

3)  Each  constrained  parameter  is  converted  to  two  con- 

constraints  in  CONMIN  unless  ABS(BL)  or  ABS(BU)  exceeds 
1.0E+15,  in  which  case  no  constraint  is  created  for  that 
bound. 
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DATA  BLOCK:   P     Omit  if  NSV=0  in  block  B 
DESCRIPTION:   SENSITIVITY  OBJECTIVES 
FORMAT:   2110 


NSOBJ  IPSENS 


NOTE:   TWO  OR  MORE  CARDS  ARE  READ  HERE 


FIELD  CONTENTS 


NSOBJ:   Number  of  separate  objective  functions 
to  be  calculated  as  function  of  the  sensi- 
tivity variables. 

IPSENS:   Print  control.   If  IPSENS. GT . 0 ,  de- 
tailed print  will  be  called  at  each  step  in 
the  sensitivity  analysis.   DEFAULT=No 
print . 


74 


DATA  BLOCK 


(CONTINUED) 


DESCRIPTION 


FORMAT 


8110 


NSN1 


NSN2 


NSN3 


NSN1+ 


FIELD 


CONTENTS 


1-i 


NSNI :   Global  variable  number  associated  with 
the  sensitivity  objective  functions. 


REMARKS 


1)   More  than  eight  sensitivity  objectives  are  allowed, 
data  cards  as  required  to  contain  data. 


Add 
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DATA  BLOCK:        Omit  if  NSV=0  in  block  B 
DESCRIPTION:   SENSITIVITY  VARIABLES 
FORMAT:   2110 


ISENS   NSENS 


NOTE:   READ  ONE  SET  OF  DATA  FOR  EACH  OF  THE  NSV  SENSITIVITY 
VARIABLES.   TWO  OR  MORE  CARDS  ARE  READ  FOR  EACH  SET 
OF  DATA. 


FIELD  CONTENTS 

1  ISENS:   Global  variable  number  associated  with 

the  sensitivity  variable. 

2  NSENS:   Number  of  values  of  this  sensitivity 

variable  to  be  read  on  the  next  card. 
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DATA  BLOCK:   Q     (CONTINUED) 
DESCRIPTION:   8F10 


FORMAT 


SNS1    SNS2    SNS3    SNS4 


FIELD  CONTENTS 

1-8  SENSI :   Values  of  the  sensitivity  variable. 

I=1,NSENS.   1=1  corresponds  to  the  nominal 
value . 


REMARKS : 

1)   More  than  eight  values  of  the  sensitivity  variable  are 
allowed.   Add  data  cards  as  required  to  contain  the  data 
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DATA  BLOCK:   V 


DESCRIPTION:   COPES  DATA  'END'  CARD 


FORMAT:   3A1 


END 


FIELD 


CONTENTS 


The  word  'END'  in  columns  1-3 


REMARKS 

1)  This  card  MUST  appear  at  the  end  of  the  COPES  data. 

2)  This  ends  the  COPES  input  data. 

3)  Data  for  the  user-supplied  routine,  ANALIZ ,  follows  this 
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HEAT  PIPE  ANALYSIS 

Data  for  the  heat  pipe  analysis  follows  the  'END'  card 
in  the  COPES  data  deck.   If • the  general  design  capability  of 
COPES/CONMIN  is  not  needed,  the  heat  pipe  analysis  can 
be  run  by  setting  NCALC  =  1  in  field  number  1  of  data  block 
B ;  or  in  a  stand-alone  mode  by  using  the  following  main 
program. 


C   MAIN  PROGRAM  FOR  HEAT  PIPE  ANALYSIS 
C   READ,  EXECUTE,  AND  PRINT  RESULTS 

DO  10  ICALC  =1,3 
10   CALL  ANALIZ  (ICALC) 

STOP 

END 
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DATA  BLOCK:   AA 


DESCRIPTION:   ELEMENT  CONNECTIVITY 


FORMAT:   315 


NEL    NSNP    NBAN 


FIELD  CONTENTS 

1  NEL:   Number  of  elements 

2  NSNP:   Number  of  system  nodal  points 

3  NBAN:   System  band  width 
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DATA  BLOCK:   BB 

DESCRIPTION:   ELEMENT  CONNECTIVITY 

FORMAT:   415  I  =  1,3;  IEL=1,  NEL 


IEL   ICOR(IEL.I) 


FIELD  CONTENTS 

1  IEL:   The  element  number 

2  ICOR(IEL,I):   System  nodal  point  corresponing 

to  nodal  point  I  of  element  IEL 


REMARKS: 

1)  Number  all  elements  with  convective  boundaries  first  from 
top  to  bottom,  then  number  the  remaining  elements. 

2)  Number  the  nodal  points  of  each  element  moving  in  a 
counter-clockwise  direction. 

3)  The  elements  with  convective  boundaries  have  nodal  points 
1  and  2  located  on  the  convective  boundary . 


DATA  BLOCK:   CC 


DESCRIPTION:   CONDENSER  GEOMETRY 


FORMAT:   7G10 . 5 


CLI    CANGL  RBASEI    R2I   THICKI   BFIN     TZ 


FIELD 


CONTENTS 


1 
2 
3 

5 
6 
7 


CLI:   Condenser  length  (inches) 
CANGL:   Cone  half  angle  (degrees) 
RBASEI:   Inside  radius  of  condenser  base 

(inches ) 
R2I:   Intermediate  radius  (inches) 
THICKI:   Condenser  wall  thickness  (inches) 
BFIN:   Height  of  fin  (inches) 
TZ :   Nodal  point  temperature  initial  guess 

(degrees  F) 


REMARKS 

])   Set  TZ  equal  to  TSS-Saturation  temperature  of  the 
working  fluid 
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DATA  BLOCK:   DD 

DESCRIPTION:   FINITE  ELEMENT  GEOMETRY 

FORMAT:   515 


NDIV    NEST    NEFB    NBOTI   NBOTF 


FIELD  CONTENTS 

1  NDIV:   Number  of  increments  along  the  length 

of  the  condenser. 

2  NEST:   Number  of  the  element  on  the  right  end 

of  the  trough . 

3  NEFB:   The  element  number  with  convective 

boundary  located  at  the  base  of  the  fin. 

M-              NBOTI:   The  element  number  with  convective 

boundary  located  at  the  right  hand  of  the 
bottom  side. 

5              NBOTF:   The  element  number  with  convective 

boundary  located  at  the  left  hand  of  the 
bottom  side. 
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DATA  BLOCK:   EE 


DESCRIPTION:   DATA  FOR  RUNNING 


FORMAT:   4F10.2 


RPM    TSS    TINF    HINF 


FIELD  CONTENTS 

1  RPM:   Rotation  rate  of  heat  pipe  (RPM) 

2  TSS :   Saturation  temperature  of  the  working 

fluid  (degrees  F) 

3  TINF:   Outside  temperature  (degrees  F) 

M-  HINF:   Outside  convective  heat  transfer 

coefficient  (BTU/HRFT  °F) 
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DATA  BLOCK:   FF 


DESCRIPTION:   CONVERGENCE  CRITERION 


FORMAT:   G10.9 


CRIT 


FIELD  CONTENTS 


CRIT:   Convergence  criterion  on  finite  element 
solution  for  incremental  heat  transfer  rate 
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DATA  BLOCK:   GG 

DESCRIPTION:   INTERNAL  FIN  GEOMETRY 


FORMAT:   2G10.5 


FANGL    ZOA 


FIELD  CONTENTS 

1  FANGL:   Fin  half  angle  (degrees) 

2  ZOA:   Ratio  of  trough  width  to  fin  base  width 
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DATA  BLOCK:   HH 


DESCRIPTION:   INTERNAL  FIN  GEOMETRY 


FORMAT:   15 


IFF 


FIELD  CONTENTS 


IFF:   (n-1) ,  where  n  is  the  number  of  rows  of 
the  upper  triangular  fin  section 


NOTE:   See  Figure  A-l 
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DATA  BLOCK:   II 

DESCRIPTION:   INTERNAL  FIN  GEOMETRY 

FORMAT :   1615  1=1, IFF 

12       3       4       5       1 


KFIN(I)  KFF(I) 


FIELD  CONTENTS 


KFIN:   The  number  of  system  nodal  points  located 
on  the  symmetric  boundary  of  triangular  fin 
section,  but  does  not  include  the  system 
nodal  points  located  at  the  base  of  the  fin 
and  the  apex. 

KFF:   The  number  of  system  nodal  points  located 
along  the  fin  convective  boundary,  but  does 
not  include  the  system  nodal  points  located 
at  the  base  of  the  fin  and  the  apex. 


NOTE:   See  Figure  A-l 
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DATA  BLOCK:   JJ 

DESCRIPTION:   INTERNAL  FIN  GEOMETRY 

FORMAT:   2G10.5,4I5 


• 


2       3       4       5       6       7 


DOBF    DOTH     JTC     JLC    JINT     KT 


FIELD  CONTENTS 

1  DOBF:   Number  of  column  within  fin. 

2  DOTH:   Number  of  column  within  though. 

3  JTC:   The  number  of  the  system  nodal  point 

located  at  the  junction  of  the  symmetry 
boundary  and  the  line  of  intersection 
between  the  fin  and  the  condenser  wall. 

4  JLC :   The  number  of  the  system  nodal  point 

located  at  the  center  of  system  coordinates . 

5  JINT:   The  numerical  difference  between  the  two 

adjacent  system  nodal  points  vertically 
at  the  condenser  section. 

6  KT:   The  number  of  rows  within  the  wall  section 


NOTE:   See  Figure  A-l. 
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IFF  =  2 
KFIN  =  2  ;  4 
KFF  =  3  ;  6 


DOBF  .'  =  3.0 
DOTH  =  i.o 
JTC   =  7 


JLG  =  17 
JINT"  =  3 
KT  =  2 


Figure  A-l.   Specification  for  Input  Data  to  Determine 
the  Coordinates  of  the  System  Nodal  Points 
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APPENDIX  B 
COMPUTER  PROGRAM  LISTING 


C  ******************************************************** 

c  *  * 

C  *   ANALYSIS  CF  ROTATING  HEAT  PIPE  ,  USIN3  TRIANGULAR    * 

C  *  ELEMENT  MODEL  * 

C  *   COMPILED  BY  MAJOR  I GNATI US. S. PURNOMO  IN  JUNE  1978    * 

C   *       MODIFIED  TO  PERMIT  NUMERICAL  OPTIMIZATION        * 

C   *  USING  COPES/CONMIN  * 

C   *     BY  LCDR  WILLIAM  A.  DAVIS,  JR.  IN  SEPTEMBER  1980    * 

C   *  * 

r   ******************************************************** 

C 

C 

C      SUBROUTINE  ANALIZ  CHANGES  THE  ARRAY  3F  DESIGN 
C      VARIABLES  FROM  SINGLE  TO  DOUBLE  PRECISION  AND 
C      BACK.   COPES/CONMIN  USES  SINGLE  PRECISION  ONLY; 
C      DOUBLE  PRECISION  IS  MAINTAINED  IN  SUBROUTINE  FUN 
C      TO  ALLOW  FOR  POSSIBLE  ILL  CONDITIONING. 
C 

c 

SUBROUTINE    ANALIZ    (ICALC) 

COMMON    /GLCBCM/    ARRAY (750) 

COMMON    /GLCB1/    BARAY(50) 

REAL*8    BARAY 

IF    (ICALC. GT.l)    GO    TO    20 

DO    10    1=1,50 

EARAY( I)=0.ODO 
10    CONTINUE 
20    CONTINUE 

CO   30    1=1,50 
30    BARAYU  ) =DBLE ( ARRAY ( I )  ) 

CALL    FUN    (ICALC) 

DO    40    1=1,50 

ARRAYd  )=SNGL(BARAY(  I)  ) 
40    CONTINUE 

RETURN 

END 
C  GUIDE    TO    FORTRAN    VARIABLE    NAMES 

C 

C         ALFA  FIN    HALF    ANGLE     (RADIANS) 

C         BFIN  HEIGHT    CF    FIN    (INCHES) 

C  BVIN  HEIGHT    OF    FIN    (FEET) 

C         CALFA         COSINE    OF    ALFA 
C         CANGL         CONE    HALF    ANGLE    (DEGREES) 

C         CBASE         INSIDE    CIRCUMFERENCE    OF    CONDENSER    (FEET) 
C         CEXIT  INSIDE    CIRCUMFERENCE    AT    CONDENSER    EXIT    (FEET) 

C         CF  THERMAL    CONDUCTIVITY    OF    CONDENSATE    FILM     (BTU/HR) 

C         CL  CONOENSER    LENGTH    (FEET) 

C         CLI  CCNCENSER    LENGTH    (INCHES) 

C         CPHI  COSINE    OF    PHI 

C         CRIT  CONVERGENCE    CRITERION 
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c 

OIV 

c 

DMTOT 

c 

DOBF 

c 

DOTH 

c 

FAMGL 

c 

H 

c 

HFG 

c 

IEL 

c 

JLC 

c 

c 

JTC 

c 

c 

c 

c 

KFF 

c 

c 

KFIN 

c 

c 

c 

KT 

c 

NBAN 

c 

NBOTF 

c 

NBOTI 

c 

NDIV 

c 

NEFB 

c 

NEL 

c 

NEST 

c 

NSNP 

c 

PHI 

c 

PI 

c 

RBASE 

c 

RBASEI 

c 

REXIT 

£ 

SALFA 

SPHI 

c 

THICK 

c 

THICKI 

c 

TPHI 

c 

ZFIN 

c 

ZGA 

c 

c 

c 

SUBR 

c 

TRAN 

c 

c 

FLOATING  POINT  VALUE  OF  ND 
CONDENSATE  MASS  FLOW  RATE 
NUMBER  OF  COLUMN  WIWTHIN  T 
NUMBER  OF  COLUMN  WITHIN  TH 
FIN  FALF  ANGLE  (DEGREES) 
CONVECTIVE  HEAT  TRANSFER  C 
LATENT  HEAT  OF  VAPOR  IZATIO 
THE  ELEMENT  NUMBER 
NUMBER  OF  SYSTEM  NODAL  POI 
THE  CENTER  OF  SYSTEM  COORD 
NUMBER  OF  SYSTEM  NODAL  POI 
THE  JUNCTION  OF  THE  SYNNET 
THE  LINE  OF  INTERSECTION  B 
AND  THE  CONDENSER  WALL 
NUMBER  OF  SYSTEM  NODAL  POI 
THE  FIN  CONVECTIVE  BOUNDAR 
NUMBER  OF  SYSTEM  NOCAL  PO 
SYSMMETRIC  BOUNDARY  OF  TR I 
NOT  CCUNTING  POINTS  AT  BAS 
NUMBER  OF  ROWS  WITHIN  THE 
SYSTEM  BAND  W IDTH 
LAST  ELEMEMT  AT  BOTTCM  SID 
FIRST  ELEMENT  AT  BOTTOM  SI 
NUMBER  OF  INCREMENT 
ELEMENT  NUMBER  AT  BASE  OF 
NUMBER  OF    ELEMENTS 
ELEMENT  NUMBER  AT  END  OF  T 
NUMBER  OF  SYSTEM  NODAL  POI 
CONE  hALF  ANGLE  (RADIANS) 
PI 

INSIDE  RADIUS  OF  CONDENSER 
INSIDE  RADIUS  OF  CCNCENSER 
INSIDE  RADIUS  OF  CONDENSER 
SINE  OF  ALFA 
SINE  CF  PHI 

CONDENSER  WALL  THICKNESS  ( 
CONDENSER  WALL  THICKNESS  ( 
TANGENT  OF  PHI 
NUMBER  OF  FINS 
RATIO  CF  TROUGH  WIDTH  TO 


IV 

HE  FIN 
E  TROUGH 

OEFFICIENT  (BTU/HR  FT2) 
N  (BTU/LBM) 

NT  LOCATED  AT 

INATE 

NT  LOCATED  AT 

RY  BOUNDARY  AND 

ETWEES  THE  FIN 

NTS  LOCATED  ALONG 

Y 

INTS  LOCATED  CN  THE 

ANGULAR  FIN  SECTION 

E  AND  APEX 

WALL  SECTION 

E 
DE 

FIN 

ROUGH 
NTS 


BASE  (FEET) 
BASE  (INCHES) 
EXIT  (FEET) 


FEET) 

INCHES) 


FIN  BASE  WIDTH 


>E  FUN  READS  INPUT  DATAt  PERFORMS  HEAT 
ANALYSIS*  AND  PRINTS  RESULTS. 
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SUBROUTINE    FUN    (ICALC) 
IMPLICIT    REAL*8(A-H,0-Z) 

COMMON  /C-LOB1/  CL I , CANGL, RBASE I, R2I , TH I CKI , BF IN, TZ , TSS 
1,TINF,HINF,FANGL,20A,ZFIN,B0A,QT0T,RPM 

01 MENS  ION  Z( 200), EPS ( 200) , HZ (200) ,XC0F(5)  ,C0F(5) ,T(200 
1),QB(200) XDMDGT(200) , U F( 200) , CF( 200) , CW( 200 ) ,AMT3T( 200 
2),R(200)  ,CINC(200)  , TB(200)  ,TT(200> , TI B (200 ) , QT INC ( 200 ) 
3,NC(200),QTOTAL(100),TE(20  0) ,R00TR(4) ,ROOTI (4) , DEL (200 
4),TBM(200)  ,BH0F(200) 

COMMON  /APOL/  DOBF ,DOTH , KF IN( 50) ,KFF ( 50) , I FF , JTC, JLC, J 
1INT,KT 

COMMON    /MAFQ/    A ( 200 ,50 )  ,F ( 200 , 1 )  , H (200 ) , TS (200  ) , TS AT, C 
1K,NEL,NSNP,NBAN,IC0R( 200,3) 
COMMCN    /PCRO/    X(200)  ,Y(200),  EZERCBVIM, THICK,  TALFA, APS 
HDEN(A1,B1,22)=(-1.0DO*(A1*22**3/3.0  00+B1*ZZ**2/2.0DO) 
I) 
C  IF    ICALC=1    READ    INPUT    CATA 

IF    (ICALC. GT.l)    GO    TO    10 
C 
C 

C  *****         INPUT    MODE  ***** 

C 

c 

C  ELEMENT    CONNECTIVITIES 

C 

READ    (5,420)    NEL,NSNP , NBAN 

WRITE    (6,430)    NEL , NSNP ,NBAN 

READ    (5,440)     ( I  EL , ( I  COR ( I  EL, I ) , I  =  1 , 3 ) , I EL=1 , N EL ) 

WRITE    (6,450) 
C  WRITE (6, 251)  ( I  EL, ( I  COR ( I  EL , I ) , 1= 1, 3 ) , I EL= 1,NEL ) 

C 

C  THE    CONDENSER    GEOMETRY 

C 

READ    (5,46C)    CL I ,C ANGL  ,RBA SEI  ,R2 I ,THI C KI , BF I N,TZ 

WRITE    (6,470)    CL I,  CANGL, RB ASE I ,R2I , THI CKI , BF I N, TZ 

READ    (5,480)     NDI V, NEST ,NEFB, NBOT I , NBOT F 

WRITE    (6,490)    NDI V, NE ST, NEFB ,NBOTI ,NBO TF 
C 

C  OATA    FOR    RUNNING 

C 

READ    (5,500)     RPM, TSS ,T INF, HINF 

WRITE    (6,510)    RPM,TSS,TINF,HINF 
C 

C  THE    CCNVERGENCE    CRITERIAN 

C 

READ    (5,520)    CRIT 

WRITE    (6,530)    CRIT 
C 

C  INTERNAL    FIN    GEOMETRY 

C 

READ    (5,540)    FANGL,ZOA 

WRITE    (6,550)    FANGL,ZOA 

READ    ( 5,560)     IFF 

WRITE    (6,570)     IFF 

READ    (5,580)     ( KFI N ( I ) ,KFF ( I ) , I =1 , I FF) 

READ    (5,590)    DOBF, DOTH, JTC  , JLC, J  INT, KT 

NHB=NEFB/2 

NTM=NBOTI+(NBOTF-NBOTI ) /2 

NBF=N30TF-H 

WRITE  (6,600)  ICOR(NB0TI ,2) ,ICOR(NEFB, 1), IC0R(NTM,2), I 
1C0R(NEST, 1),IC0R(NB0TF,1) 

RETURN 
C  IF    ICALC=2    PERFORM    THE    HEAT    TRANSFER    ANALYSIS 

10     IF    (ICALC. GT. 2)    GO    TO    360 
C 
C 
C  *****         EXECUTION    MODE         ***** 
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c 
c 

C  CONVERT    UMTS    OF    ALL    DIMENSIONAL    PARAMETERS 

C  TO    FEET.       CONVERT    UNITS    OF    ANGLES    TO    RADIANS. 

CL=CLI/12.0D0 

R2=R 21/12. ODO 

RBASE=RBASEI/12.0D0 

BVIN=BFIN/12.0D0 

01  V=DFLOAT(NOIV) 

P 1=3. 141 5926535 8979 DO 

PHI=2.0DO*CANGL*PI/3  60.0DO 

SPHI=DSIN(PHI ) 

CPHI=DCOS(PHI) 

TPHI=0TAN<PHI ) 

CELX=CL/CIV 

CBASE=2.0D0*PI*RBASE 

REXIT=RBASE+CL*TPHI 

CEXIT=2.0D0*PI*REXIT 

THICK=THICKI/12.0D0 

ALFA=FANGL*2.0D0*PI/360. 00 0 

SALFA=DSIMALFA) 

CALFA=OCOS(ALFA) 

TALFA=DTAN(ALFA) 

EZERO=2.0CO*BVIN*TALFA 
C 

C  BOUNCARY    CONDITIONS    AND    TEMPERATURE    ESTIMATES 

C  ALCNG    THE    FIN    BOUNDARY 


C 


DO    20    NTINF=NBOTI, NBOTF 
20    TS(NTINF)=TINF 

00    30    NNT=NBF,NEL 

TS(NNT)=0.0D0 
30    H(NNT)=C.0DC 

DO    40     IGT=1,NEST 

IE  =  IC0RUGT,2) 
40    T(  IE)=TZ 

IG=IC0R(NEST,1) 

T(IG)=TZ 

CMEGA=RPM*2.0D0*P 1*60.  000 

DO    50    KL=NBCTI,NB0TF 
50    H(KL)=HINF 

FIFN=HINF 

TSAT=TSS 

EPSO=ZGA*EZERO 

E0A=BVIN/(EZER0/2.0D0) 

ZFIN=CBASE/(EZERO+EPSC) 

SURFAR=ZFIN*(2.0D0*( BV  IN /CALF A ) +EP SO ) 

EPSEX=<CEXIT-(ZFIN*EZERO) )/ZFIN 

BETA=(  EPSEX-EPSOJ/DI V 

ZZERG=BVIN/CALFA 

ZA=0.000 

DO    60    NSAT=1,NEST 
60    TS(NSAT)=TSAT 

TSQLID=( TSAT+TINF) /2.0D0 
TEMPORARY    CHANGE    -    TFILM 

GT=O.ODO 

QBTOT=0.CDC 

QT1=0.000 

QTF^C.ODO 

QTOT=0.000 

DMTOT=0.0D0 

NK=NDI V+l 

DO    350    NI=1,NK 

R(NIi  =  R2  +  M*DELX*SPHI 

EPS(NI )=EPSG+NI*BETA 

APS=EPS(NI  ) 
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C  NCCAL    POINT    COORDINATES 

CALL    COORD 

Z(1)=ZA 

00    70    IZEL=1,NEFB 

NA=IC0R(IZEL,  1) 

NB=IC0R(IZEL,2) 

XE=X(NA)-X(NB) 

YE=Y(NA)-Y(NB) 

ELZ=DSQRT(XE**2+YE**2) 
70    U  IZEL+1)=Z( IZELJ+ELZ 

XZB=XUCCP(NHB,1)  )-X(  ICOR(  1,  2)) 

YZB=Y(ICCRCNHBtl)  ) -Y ( I  COR < 1, 2 ) ) 

ZB=OSQRT (XZB**2+YZB**2) 

ZC=ZZERO 

IM=1 
C 

C  PARABOLIC    TEMPERATURE    DISTRIBUTION    ALONG   THE    FIN 

C  BOUNDARY, USING    LAGRANGE    INTERPOLATION 

C 

80    TP1  =  TUCCR(1,2)  ) 

TP2=T<  ICOR(NHBt 1) ) 

TP3=T(ICCR<NEFB,1 ) ) 

AP1=TP1/(ZB*ZC) 

AP2=TP2/ (ZB*(ZB-ZC) ) 

AP3=TP3/(ZC*(ZC-ZB> ) 

BP1=-(ZB+ZC)*AP1 

BP2=-ZC*AP2 

BP3=-ZB*AP3 

A1=AP1+AP2+AP3 

B1=BP1+6F2+BP3 

TC=0.000 

CO    90    NY=1,NEST 
90    TC=TC+T(ICCR(NY,2) ) 

AY=DFL0AT(NY+1) 

TF=(TC+T(ICCR(NY, 1) )+AY*TS(NY) )/(2.0D0*AY) 
C 

C  SOLID-FLUID    PROPERTIES 

C 

HFG=109  7.2DO-0.60187  50  0*TS(1) 

RHOF (N I) =62.77400-0. 00 255698D0*TF- 0.00 005 357 2D0*TF** 2 

CF(NI) =0.3 03400+0. 0  00 73892 7D0*TF-0. 000 00 14732 1D0*TF**2 

UF(NIJ=0.0013  97D0-0. 00 00 146690 0*TF+ 0.00000006 31 25 3D0*T 
lF**2-0 .00000000009 76 56 9D0*TF**3 

LF(NI) =36CC*UF(NI ) 
C      CW ( NI )  =  23 1.777200- 0.02 22 2D0*T SOL  ID 

CW(NI)=8.776+0.0026  5D0*TSOLID 
C      CW(NI)=1.C 
C      CW(NIJ  =  20000.0 

CK=CW(NI) 

CONST=RHOF(NI)**2*OMEGA**2*HFG*CPHI*CALFA*R(NI) 
C 

C      INITIAL  FILP  THICKNESS 
C      DEL(1)=0.00006752DO 

IF  (NI.GT.l)  GO  TO  100 

OEL( 1)=1.1C7*(( (TSAT-TINF)*CF ( NI ) / { UF ( NI ) *HFG*3600 .0  00 
1 ))**0.25)*(<UF(NI)/(RH0F(NI)*0MEGA) )**0.5) 
100  CONTINUE 
C 

C  AVERAGE  ELEMENT  CONVECTIVE  COEFFICIENT  ALONG 

C  THE  FIN  BOUNDARY 


C 


ZSTAR=ZZERC-DEL(NI  )/CALFA 

AZZ=DEL(NI)/SALFA 

ZZ=ZSTAR 

AZS=DABS(4*CF(NI)*UF(NI)*HDEN(A1 , Bl ,  ZZ )/CCNST )**0.25D0 

HAC=O.ODO 
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00    190    IEL=1,NEFB 

AZ=Z(IEL) 

BZ=Z(IEL>1) 

IF     (ZSTAR.LE.BZ)    GO    TO    110 

GO    TO    120 
110    IF    (HAC.NE.O.ODO)     GO    TO    180 

8Z=ZSTAR 
120     IF    (  IEL.NE.1)    GO   TO     130 

AK={3Z-AZ)/5.0D0 

ZZ=AK 

GO    TO    140 
130    AK=(6Z-AZ)/4.0D0 

2Z=AZ 
140    ZEL=4*AK 

00    150    Nh=l,5 

HZ(NH)=DABS(CF<NI)**3*C0NST/(4*UF(NI)*HDEN(A1,B1,ZZ) ) ) 
1**0.2500 
150    ZZ=ZZ+AK 

CGNH=AK*(HZU)+4*HZ(2)+2*HZ(3)+4*HZ(4)  +HZ  (  5  )  )  /  {  3*Z  EL) 

IF    (ZSTAR.EQ.BZ)    GO    TO    160 

H(IEL)=CCNH 

GO  TC  190 
160  AZ=ZSTAR 

HAZ=CONH*<AZ-Z( I EL ) ) 

DELA=AZS 
170    BZ=Z(IEL+1) 

DELB=(8Z-ZST;*R)*AZZ/ (ZZERO-ZSTAR ) 

0tLZ=(DELA+D£LB)/2.0D0 

FAC= (BZ-AZ)*CF(NI)/DELZ 

HUEL)=(HAZ+HAC)/(BZ-Z(IEL)) 

GO    TO    190 
180    AZ=Z(IEL) 

DELA=OELB 

HAZ=O.ODO 

GO    TO    170 
190    CONTINUE 

NETI=NEFB+1 

DO    200    IEL=NETI,NEST 
200    M  IEU=CF(NI)/DEL(NI) 
C 

C  ENTRY    INTO    THE    FINITE    ELEMENT    SOLUTION 

C 

CALL    FCRMAF 

CALL    BANOEC     ( NSNP, NB AN , 1 ) 
C 

C  THE    TEMPERATURE    DISTRIBUTION 

C 

CO    210    NT=1,NSNP 
210    T(NT)=F(NT,1) 

TIB(N1)=T( IC0R(NB0TI,2)) 

TT(NI)=T(ICCR(NEFB, 1) ) 

TBM(NI )=T(IC0R(NTM,2) ) 

TE(NI)=T(ICOR(NEST, 1)  ) 

TB(NI)*T(ICCR(NB0TF,1 ) ) 

TTS=O.OOC 

DO    220    N5=1,NSNP 
220    TTS=TTS+T(NS) 

PN=DFLOAT(NS) 

TSOLID=TTS/PN 
C 

C  Q    AT    THE    BOTTOM    SIDE 

C 

QBI=G.OOC 

DO    230     IEEL=NBOTI,MBOTF 

NKA=IC0R(IBEL,1) 

NKB=IC0R(IBEL,2) 

XB=X(NKA)-X(NKB) 
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YB=Y(NKA)-Y(NKB) 

EL B-DS QRT ( X B**2 +Y B**2 ) 
230    GBI=QBI+<T(NKA)+T(NKB)-2*TS( IB  EL ) ) *ELB*H ( I  BEL ) /2.0 00 

CB(NI)=QBI*DELX 
C 

C  ITERATION    UNTIL    CONVERGENCE    CRITERIA    IS    MET 

C 

IF    { IM.EQ.l)    GO   TO    240 

QJ=QBI 

GO   TO    250 
240    QI=QBI 

IM=2 

GO    TO    80 
250    AQ=DABS(CJ-CI)/QJ 

IF    (AG.LE.CRIT)    GO    TO    260 

CI=QJ 

GO    TO    80 
260    OMDOT(NI )=2.*QBI*DELX/HFG 

CMTQT=DMTCT+DMDOT(NI ) 

Cl=RHCF(M)**2*OMEGA**2*R(NI)*SPHI/(3*UF(NI)  ) 

XCOF(l)=-Df*TOT 

XCOF(2)=0.CC0 

XCOF(3)=0.0DO 

XC0F(41=C1*EPS(NI) 

XC0F(5  )=C1*TALFA 

M=4 

CALL    DPOLRT    ( XCOF , COF , M, ROCTR, ROOTI , I ER) 

IF    (RGOTP(l) .GT.O.ODO)    GO    TO    270 

IF    (R00TR(2).GT. 0.000)    GO    TO    280 

IF    (RG0TR(3)  .GT.O.ODO)    GO    TO    290 

IF    (RG0Tf!(4)  .GT.O.ODO)    GO    TO    300 

kRITE    (6,610) 

WRITE    (6,620)     ( ROOTR (  I  ) , 1= 1, 4 ) 
C  GO    TO    6100 

C 

C  THE    CONDENSATE    THICKNESS 

C 

270    CEL(NI+1)=R00TR(1) 

GO    TG    310 
280    DEL(NH-1)=R00TR(2) 

GO    TO    310 
290    0EL(NI+1)=RCGTR(3) 

GO    TO    310 
300    0EL(NI-H)  =  R00TR(4) 
310    OEL=O.OOG 

IF    (NI.NE.l)    GO    TO    320 
C 

C  Q    FROM    THE    TCP    SIDE 

C 

C  Q    THROUGH    FIN 

C 

320    DO    330    ICEL=1,3 

KA=ICOR(ICEl,l) 

KB=IC0R(ICEL,2) 

XQEL=X(K4)-X(KB) 

YQEL=Y(KB)-Y(KA) 

ELM=DSQRT(XCEL**2+YQEL**2) 

QEL=QEL+(2*TS(IQEL)-T{KA)-T(KB) ) *ELM*H ( IQEL ) /2 .000 
330    CONTINUE 

CINC(NI)=QEL*DELX 

AMTOT(NI )=OMTOT 

CET=QEL*CELX*ZFIN*2 

QT=QT+QET 

QA=QBI*DELX*ZFIN*2 

QTOT=QTOT+QA 
C 
C  C    THROUGH    TROUGH 


97 


DO    340    IQEL=4,4 

KA=IC0R(ICEL,1) 

KB=IC0R(IQEL,2) 

XQEL=X (KA)-X(KB) 

YQEL=Y(K8)-Y(KA) 

ELM=DSQRT( >QEL**2+YQE L**2> 

CTRF=(2*TS(  IQEl)-T(KA)-T(KB)  )*ELM*H(  I3EU/2.0D0 
340    CONTINUE 

CTINC(NI)=QTRF*DELX 

CTGT AL ( NI ) =CINC(NI)+QT IN C(NI) 

QTRFT=QTRF*0ELX*ZFIN*2. 

CTF=QTF+QTRFT 
350    CONTINUE 

RETURN 
360    CONTINUE 
C 

c 

C  *****         CUTPUT    MODE         ***** 

WRITE     (6,630) 
CO    370    NR=1,NDIV 
370    WRITE    (6,640)    NR, QI NC ( NR ) , GT INC( NR ) , QT OT AL (NR ) 
WRITE    (6,650)    QT,QTF 

WRITE    (6,660)    CL I, CANGL, RBASE I,R 21 , THI CK I , BF IN ,RPM , TSS 
1,TINF,HINF,CRIT,FANGL,Z0A,  IFF 
WRITE    (6,670)    BOA , ZOA , ZFI N ,BVIN, SURFAR 
WRITE    (6,680) 
DO    380    NP=1,NSNP 
380    WRITE    (6,690)    NP, X ( NP ) , Y( NP ) , T(NP ) 
WRITE    (6,700) 
CO    390    KKL=1,NB0TF 
NKX=IC0R(KKL,1) 
NKY=IC0R(KKL,2) 
XP=X(NKXJ-X(NKY) 
YP=Y(NKX)-Y(NKY) 
EXY=0SQRT(XP**2*YP**2) 

CEP=DABS((T(NKX)+T(NKY)-2*TS(KKL) )*EXY*H(KKL) /2.0D0) 
CEP=QEP*CELX 
390    WRITE    (6,710)    KKL,H ( KKL) ,EXY , QEP 
WRITE    (6,720)    CRIT 

WRITE     (6,730)    HFG , ZFI N ,H (NBOTF ) , TS AT ,R PM, QTOT , QT, FANGL 
WRITE    (6,740) 
00    400    NR=1,NDIV 
400    WRITE    (6,75C)    NR,DEL ( NR) , QE( NR) , AMTOT ( NR ) ,T I B ( NR) , TT ( N 
1R),TE(NR),T8(NR) 
WRITE     (6,760) 
DO    410    NG=1,NDIV,2 
410    WRITE    (6,770)    NG, CW ( NG  ) , CF (NG ) , RHOF ( MS ) , UF ( NG ) ,EP S( NG ) 
1,R(NG) ,TBW(NG) ,QINC(NG) 
RETURN 
C 

420    FORMAT    (315) 

430  FORMAT  (/2X, 15HN0 .0 F . ELEMENTS=, I  5 , 10X , 34HNC. OF. SYSTEM 

1  N.P.=,I5,10X,13HN0.0F  BANDED=,I5) 

440  FORMAT  (415) 

450  FORMAT  </2X, 7HELEMENT , 10X, 3HNP 1, 14X,  3HNP 2, 15X, 3HNP3) 
460  FORMAT  (7G10.5) 

470  FORMAT  (4X,5HCLI=  , E 12. 5, / ,4X, 7HCANGL=  , E12. 5  ,/ ,4X ,8HR 
ieASEI.=,E12.5,/,4X,5HR2I=  , E12.5, / , 4X, 8HTHICK 1=  ,E12.5 
2,/,4X,6HBFIN=  , E12. 5 , / ,4X,4HTZ=  ,E12.5) 
480  FORMAT  (515) 
490  FORMAT  (4X,6HNDIV=  , I 10 ,/ ,4X ,6HNEST=  , I 10,/ , 4X, 6HM EFB= 

1  ,I10,/,4X,7HNB0TI=  , I  10, / ,4X ,7HNB0TF=  ,110) 
500  FORMAT  (4F10.2) 
510  FORMAT  (4X,5HRPM=  , E12 .5  ,/ ,4X  ,  5HTSS=  , E12 .5 , / ,4X, 6HT IN 

1F=  ,E12.5,/,4X,6HHINF=  ,E12.5) 
520  FORMAT  (G10.9) 
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530  FORMAT  (4X,6HCRIT=  ,E12.5) 

540  FORMAT  (2G10.5) 

550  FORMAT  (4X,7HFANGL=  ,  612. 5 ,/ ,4X, 5HZ0A  =  ,E12.5) 

560  FORMAT  (15) 

570  FORMAT  (4X,5HIFF=  ,110) 

580  FORMAT  (  1615) 

590  FORMAT  <2G10.5,4I5) 

600    FORMAT     { ///5X  ,4HTI B  =  , I  5 ,10X,3HTT  =  , 15 , / , 5X,4HT BM=, I  5 , 10 

1X,3HTE=, I5,/,6X,3HTB=,I5) 
610    FORMAT     (// 10X , 17HCRASH, CRASH, CRASH) 
620    FORMAT    ( //5X , 4( E12. 7 , 3X) ) 
630    FORMAT     { 2X, 7hELEMENT, 2X,4HQF IN, 1 7X, 7HQ TROUGH , 15X, 6HQT0 

1TAL) 
640    FORMAT    ( 4X, I 5,E 12. 5, 1 0X,E12. 5 , 10X, E12. 5) 
650    FORMAT     (// ,4X, 11HQF IN    TOTAL= , E 12 . 5, 10X, 15HQTRCUGH    TOTA 

1L=    ,E12.5) 
660    FORMAT    (///// ,4X, 5HCL I =    , E  12. 5, 5X, 7HCANGL=    ,E12.5,/,4X 
1,8HR8ASEI=    ,E12.5,2X,5HR2I=    , E12 .5, /, 4X, 8HTHICKI=     ,E12 
2.5,2X,6HBFIN=    , E12 . 5 , / ,4X, 5HRPM=    , El  2 . 5 , 5X ,5HTSS=    , E12 
3.5,/,4X,6hTINF=    , E 12 . 5, 4X, 6HH INF=    ,E 12. 5  ,  / ,4X  ,6HCRI T= 
4,E12.5,4X,7HFANGL=    , E12 .5 , / ,4X ,5 H20A=     , E12.5, 5X, 5H I FF= 
5    ,110) 
670    FORMAT     ( 1H1 , //2X.4HB0 A=, Gl 2. 5, 5X, 4HZ0 A=, G 12. 5 , 5X, 5HZF I 

1N=,G12.5,5X,5HBFIN=,G12.5,5X,13HSURFACE    AREA=,G12.5) 
680    FORMAT     < // 5X , 2HNP, 6X, 1HX, 12X, 1HY, 12X , 1HT) 
690    FORMAT     i/ZX ,  13 , 3X, 3 ( F10 .6, 3X ) ) 

700    FORMAT    ( //2X,2HEL , 8X , 1HH, 1 IX ,9HEL-LENGTH , 15X , 4HQ-EL ) 
710    FORMAT     (  /  2X,  12,  3X,  E12  .  5,  3X  ,E  12.  5  ,  10X  ,  E  12.  5) 
720    FORMAT     ( /2X,22HCCNV ERGENCE    CR IT ER I AN= , E15  .8 ) 
730    FORMAT    ( 1H    , // , 5X, 4HHFG= ,E  12. 5 ,/ , 5X , 1 1HN0.CF    FINS=,E12 
1.5,/,5X,6FH-0UT=,E12.5,/,5X,5HTSAT=,E12.5,/, 5X,4HRPM=, 
2E12.5,/,5X,6HQ-B0T=,E12.5,/,5X,6HQFIN    =,E12.5,/,5X,11H 
3HALF-ANGLE=,F8.3) 
740    FORMAT     ( 1H0 , 6X , 1HJ, 4X, 14HF ILM    TH ICKNESS, 6X, 8HQ-INCREM , 

16X,8HMASS-TCT,7X,3HTIB,8X,2HTT,10X,2HTE,8X,2HTB) 
750    FORMAT     ( 1H    , 4X, 14,  4X, F 12. 10, 4X,F 10.4, 6X ,F 9.  5 , 6X ,F 5.  1 ,6 

1X,F5.1 ,6X,F5.1,6X,F5.1 ) 
760    FORMAT    ( 1HC, 6X, 1H J , 6X, 6HK- WALL ,4X,6HK-FI LM,3X ,7H0ENSIT 
1Y,4X,9HVISC-FILM,6X,7HEPSIL0N,5X,6HRADIUS,5X, 3HTBM,5X, 
25HQ-BCT) 
770    FORMAT     ( lh    , 4X, 14, 4X, F7. 3, AX,F6. 4, 4X ,F 6. 3 ,4X , F9.7 ,4X,F 
19.7,4X,F7.5,5X,F5.1, IX, 1P1C12.3) 
END 

SUBROUTINE    COORD 
IMPLICIT    REAL*8<A-H,0-Z) 

COMMON    /GL0B1/    CLI ,CANGL , RBA SEI , R2I , TH I CKI ,BF I N,T Z, TSS 
l,TINF,HINF,FANGL,Z0A,ZFIN,BOA,QTOT,RPM 
COMMON    /PCRO/    X(200)  ,  Y  (200  )  ,  EZERO,  BV/IN,  THICK,  TALFA,  APS 
COMMON    /APQL/    DOBF , DOTH, KF IN( 50)  ,KFF ( 50)  ,  IFF  ,  JTC, JLC , J 
1INT,KT 
DtLH=BVIN/DCBF 
X(1)=0.0C0 
Y(1)=THICK+BVIN 
N=l 

DO    20     1=1,  IFF 
ICA=KFIN(I ) 
ICB=KFF( I) 
CBA=DFLOAT(ICB-ICA) 
AN =0.000 

DO    10    II=ICAfICB 
X(II)=X(1)+N*AN*DELH*TALFA/CBA 
Y(II)=Y( 1)-N*DELH 
10    AN=AN+1.0C0 
20    N=N+1 

AN=O.ODO 

ICD=ICB-ICA+1 

DO    50    J=JTC,JLC,JINT 
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X(  J)=X(1) 

Y(  J)=( 1. COO- AN/DOTH)* THICK 

CO    30    JJ=1,  ICO 

X< J+JJ)=X(J)+JJ*EZER0/ (2*(CBA+1.0D0) J 
30    Y(  J*JJ)=:Y(  J) 

CO    4-0    K=1,KT 

X( J+JJ+K)=X< J+JJ)+K*APS/(2.0D0*KT) 
40    Y( J+JJ+K)=Yl J) 
50    ANaAN+i.OCO 

RETURN 

END 

SUBRCUTINE  FCRMAF 

IMPLICIT  REAL*8(A-H,0-Z) 

CIMENSION  E(3),C<3). EA(3,3) 

COMMON  /GLGB1/  CLI t CANGL , RBAS EI , R2 I ,TH ICK I, BF  IN, TZ, TSS 
1,TINF,HINF,FANGL,Z0A,ZFIN,B0A,QTCT,RPM 

COMMCN  /PCRC/  X(200) , Y ( 200  ), EZERO, BVIN, THICK, TALFA, APS 

COMMON  /MAFO/  A( 200, 50 ) , F (200, L )  ,H(200 ) , TS (200  )  ,TS AT , C 
1K,NEL,NSNP,NBAN, ICOR( 200,3) 

00  20  N=1,NSNP 

F(N, 1)=C.CDC 

DO    10    MA=1,N8AN 
10    A(N,MAJ=Q.CD0 
20    CONTINUE 

CO    70    IEL=1,NEL 

IA=ICOR( IEL,1) 

IB=IC0RUEL,2) 

IC=IC0R(IEL,3) 

B(  1)=Y(IB)-Y(IC) 

fi(2)=Y(IC)-Y(IA) 

B(3)=Y(IA)-Y(IB) 

C(U=X(IC)-X(IB) 

C(2)=X(IAJ-X(IC) 

C(3)=X(IB)-X(IA) 

EL=DSQRT ( C ( 3 ) **2+B  (  3 ) **2 ) 

AS=DABS(  <6(1)*C(2)-B(2)*C(1) J/2.0D0) 

HC-HUEU/CK 

DO   60    J=l,3 

JJ=ICGR(IEL,J) 

00    50    K=l,3 

KK=IC0R(IEL,K) 

EA{J,K)=(B( J)*B(K)+C (J)*C(K) )/(4*AS) 

IF     (HC.EG.O.ODO)    GO    TO    40 

HEL=HC*EL/6.0D0 

IF    (J.EQ.3)    GO    TO    40 

IF    (K.EQ.3)     GO    TO    ^0 

IF    (J.EQ.K)    GO    TO    30 

EA(J,K)=EA( J,K)+HEL 

GO   TO    40 
30    EA(J,K)=EA( J,K)*2*HEL 
40     IF    (KK.LT.JJ)    GO    TO    50 

NW=KK-JJ+1 

A(JJ,NW)=A< JJ,NW)+EA( J,K) 
50    CONTINUE 
60    CONTINUE 

FE=HC*TS<IED*EL/2.0D0 

F(  IA,1)=F(  IA.D+FE 

F(IB,1)=F(IB,1)+FE 
70    CONTINUE 

RETURN 

END 

SUBRCUTINE    BANDEC    (N EQ ,MAX3, NVEC ) 

IMPLICIT    REAL*8(A-H,0-Z) 

COMMON  /GL0B1/  CLI , C ANGL , RBA SE I , R2 I , TH ICKI ,BF IN,TZ , TSS 
1,TINF,HINF,FANGL,Z0A, ZFIN, BGA,QTGT,RP^ 

COMMON  /PCRO/  X(200) ,Y(200) , EZERO, BVIN, THICK, TALFA, APS 

COMMON  /MAFO/  A( 200, 50 ) , F( 200, 1) ,H( 200 ) , TS( 20C) ,TSAT,C 
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1K,NEL,NSNP,NBAN,ICOR(200,3> 

LOOP-NEQ-1 

DO   20    I=1,LCCP 

MB=I+1 

NB=MIN0<  I*MAXB-1,NEQ  ) 

00    20    J=MB,NB 

L=J+2-MB 

D=A( I,L)/A(I,1) 

00    10    MM=1,NVEC 
10    F(J,MM)=F(J,MM)-D*F( I , MM ) 

MM=MIN0(MAX6-L+1»NEQ-J+1) 

CO    20    K-1,MM 

NN=L+K-1 
20    A( JtKJ*A(JtK)-D*A(I,NN) 

00    30    I-l.NVEC 
30    F(NECtI)=F(NEQ, I)/A(NEQf 1) 

00    50-  I=2,NEQ 

J=NEQ-I+1 

K=MIN0(NEC-J*1»MAXB) 

CO    50    MM-1,NVEC 

DO   40    L=2,K 

MB=J+L-1 
40     F( J, MM J=F{ J,MM)-A( J, L) *F(MB,MM) 
50    F( J,MM)=F(J,MM)/A( J, 1) 

RETURN 

END 
C  SUBROUTINE    OPOLRT    COMPUTES    THE    ROOTS    OF    A    REAL 

C  POLYNOMIAL    USING   A    NE  WTON-RAPHSON    ITERATIVE 

C  TECHNIQUE. 

C 
C 

SUBROUTINE    DPOLRT     < X COF, CCF, M, ROOTR, RDOTI , IER ) 

IMPLICIT    REAL*8    (A-H),    REAL*8     (0-Z) 

DIMENSION    XCOF(  1),C0F(  1)  tRCOTR(  1)  ,ROOTI(  1) 
C 

c 

IFIT=0 

N  =  M 

IER=0 

IF    <XC0F(N+1) )     10t40, 10 
10    IF    (N)     20,20,60 
C 
C  SET    ERRCR    CODE    TO    1 


C 


20     IER=1 
30    RETURN 


C 

C  SET    ERRCR    CODE    TO    4 

C 

40     IER=4 

GO    TO    30 
C 

C  SET    ERRCR    CODE   TO    2 

C 

50    IER=2 

GO    TO    30 
60    IF    (N-36)     70,70t50 
70    NX=N 

NXX=N+1 
N2=l 
KJ1=N+1 
DO    80    L=1,KJ1 
MT=KJ1-L+1 
80    COF( MT)=XCCF(L> 
C 

C  SET     INITIAL    VALUES 

C 
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90    XO=. 00500101 
Y0=0. 01000101 

c 

C  ZERO    INITIAL    VALUE    COUNTER 

C 

IN=0 
100    X=XO 

c 

C  INCREMENT    INITIAL    VALUES   AND    COUNTER 

C 

X0=-10.0*YO 

Y0=-10.0*X 
C 

C  SET    X   ANO    Y    TO    CURRENT    VALUE 

C 

x=xo 

Y=YO 

IN=IN+1 

GO    TO    120 
110    IFIT=1 

XPR=X 

YPR=Y 
C 

C  EVALUATE    POLYNOMIAL    ANO    DERIVATIVES 

C 

120    ICT=0 
130    UX=0.0 

UY=0.0 

V=0.0 

YT=0.0 

XT=1.0 

U=CQF<N+1J 

IF    (U)    140,270,140 
140    DO    150    1=1, N 

L=N-I+1 

XT2=X*XT-Y*YT 

YT2=X*YT+Y*XT 

U=U+CCF(L)*XT2 

V=V+CCF<L)*YT2 

FI=I 

UX=UX+FI*XT*COF(L ) 

UY=UY-FI*YT*CCF{L) 

XT=XT2 
150    YT=YT2 

SUMSQ=UX*LX+UY*UY 

IF    (SUMSQ)     160,230,  160 
160    DX=( V*UY-U*UX)/SUMSG 

X=X+DX 

DY=-(U*UY*V*UX )/SUMSQ 

Y=Y+DY 

IF    (DABS(DY)+DABS(DX)-1.0E-05)     210,170,170 
C 

C  STEP    ITERATION    COUNTER 

C 

170    ICT=ICT+1 

IF    (ICT-5GC)     130,180,180 
180     IF    (IFIT)    210,190,210 
190    IF    (IN-5i    100,200,200 

C  SET     EPPCR    CODE    TO    3 

C 

200     IER=3 

GO  TO  30 
210  00  220  L=1,NXX 
MT=KJ1-L+1 
TEMP=XCOF(MT) 
XCOF(MT)=CCF(LJ 
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220    COF(L)=TEPP 

ITEMP=N 

N=NX 

NX=ITEMP 

IF  (IFIT)  250,110,250 
230  IF  (IFIT)  240,100,240 
240    X=XPR 

Y=YPR 
250    IFIT=0 

IF    (DABS(Y/X)-1.0E-04)     280,260,260 
260    ALPHA=X+X 

SUMSQ=X*X+Y*Y 

N=N-2 

GO    TO    290 
270    X=0.0 

NX=NX-1 

NXX=NXX-1 
280    Y=O.C 

SUMSQ=0.0 

ALPHA=X 

N=N-1 
290    C0F(2)=CCF(2)-»-ALPHA*C0F(l) 

00    300    L  =  2,N 
300    COF(L+l)=CCF(L+l)+ALPHA*COF(l)-SUMSa*COF(L-l) 
310    R00TI(N2)=Y 

R00TR(N2)=X 

N2=N2  +  1 

IF    (SUMSG)     320,330,320 
320    Y=-Y 

SUMS  0=0.0 

GO    TO    310 
330     IF    (N)    30,30,90 

END 
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